Translate this page:
You can translate the content of this page by selecting a language in the select box.
Disclaimer: The purpose of the Open Case Studies project is to demonstrate the use of various data science methods, tools, and software in the context of messy, real-world data. A given case study does not cover all aspects of the research process, is not claiming to be the most appropriate way to analyze a given dataset, and should not be used in the context of making policy decisions without external consultation from scientific experts.
This work is licensed under the Creative Commons Attribution-NonCommercial 3.0 (CC BY-NC 3.0) United States License.
To cite this case study please use:
Wright, Carrie and Ontiveros, Michael and Jager, Leah and Taub, Margaret and Hicks, Stephanie. (2020). https://github.com/opencasestudies/ocs-bp-co2-emissions. Exploring CO2 emissions across time (Version v1.0.0).
This case study explores how different countries have contributed to Carbon Dioxide (CO2) emissions over time and how CO2 emission rates may relate to increasing global temperatures and increased rates of natural disasters and storms. We used this report from the EPA as the basis for motivating this case study, as it provides background information about how CO2 emissions and other greenhouse gases have influenced the climate and weather patterns.
CO2 makes up the largest proportion of greenhouse gas emissions in the United States:
So why should we pay attention to greenhouse gases?
According to the US Environmental Protection Agency (EPA) Inventory of U.S. Greenhouse Gas Emissions and Sinks 2020 Report:
Greenhouse gases absorb infrared radiation, thereby trapping heat in the atmosphere and making the planet warmer. The most important greenhouse gases directly emitted by humans include carbon dioxide (CO2), methane (CH4), nitrous oxide (N2O), and several fluorine-containing halogenated substances. Although CO2, CH4, and N2O occur naturally in the atmosphere, human activities have changed their atmospheric concentrations. From the pre- industrial era (i.e., ending about 1750) to 2018, concentrations of these greenhouse gases have increased globally by 46, 165, and 23 percent, respectively (IPCC 2013; NOAA/ESRL 2019a, 2019b, 2019c).
* IPCC stands for the Intergovernmental Panel on Climate Change
In fact, there are many signs that our planet is experiencing warmer temperatures:
The connection between greenhouse gas levels and global temperatures and the influence of increased global temperatures on human health are motivated by these reports:
Melillo, J.M., T.C. Richmond, and G.W. Yohe (eds.). 2014. Climate change impacts in the United States: The third National Climate Assessment. U.S. Global Change Research Program.
The National Climate Assessment Report states that:
Heat-trapping gases already in the atmosphere have committed us to a hotter future with more climate-related impacts over the next few decades. The magnitude of climate change beyond the next few decades depends primarily on the amount of heat-trapping gases that human activities emit globally, now and in the future.
See the following links for more information about how greenhouse gases have influenced global temperatures: 1) The EPA report on green house gases
2) The National Climate Assessment (NCA) summary from 2014) 3) The World101 website about how countries are adapting to climate change
Our main questions:
In this case study, we will explore CO2 emission data from around the world. We will also focus on the US specifically to evaluate patterns of temperatures and natural disaster activity.
This case study will particularly focus on how to use different datasets that span different ranges of time, as well as how to create visualizations of patterns over time. We will especially focus on using packages and functions from the tidyverse, such as dplyr, tidyr, and ggplot2.
The tidyverse is a library of packages created by RStudio. While some students may be familiar with previous R programming packages, these packages make data science in R especially legible and intuitive.
The skills, methods, and concepts that students will be familiar with by the end of this case study are:
Data Science Learning Objectives:
dplyr for data wranglingdplyrggplot2ggplot2 plotsggplot2 plotsStatistical Learning Objectives:
We will begin by loading the packages that we will need:
library(here)
library(readxl)
library(readr)
library(dplyr)
library(magrittr)
library(stringr)
library(purrr)
library(tidyr)
library(forcats)
library(ggplot2)
library(directlabels)
library(ggrepel)
library(broom)
library(patchwork)Packages used in this case study:
| Package | Use in this case study |
|---|---|
here |
to easily load and save data |
readxl |
to import the excel file data |
readr |
to import the csv file data |
dplyr |
o view and wrangle the data, by modifying variables, renaming variables, selecting variables, creating variables, and arranging values within a variable |
magrittr |
to use and reassign data objects using the %<>%pipe operator |
stringr |
to select only the first 4 characters of date data |
purrr |
to apply a function on a list of tibbles (tibbles are the tidyverse version of a data frame) |
tidyr |
to drop rows with NA values from a tibble |
forcats |
to reorder the levels of a factor |
ggplot2 |
to make visualizations |
directlabels |
to add labels to plots easily |
ggrepel |
to add labels that don’t overlap to plots |
broom |
to make the output form statistical tests easier to work with |
patchwork |
to combine plots |
The first time we use a function, we will use the :: to indicate which package we are using. Unless we have overlapping function names, this is not necessary, but we will include it here to be informative about where the functions we will use come from.
Now we will describe a bit more background about greenhouse gas emissions and the potential influence of these emissions on public health.
Greenhouse gas emissions are due to both natural processes and anthropogenic (human-derived) activities.
These emissions are one of the contributing factors to rising global temperatures, which can have a great influence on public health as illustrated in the following image:
According to the US Environmental Protection Agency (EPA) Inventory of U.S. Greenhouse Gas Emissions and Sinks 2020 Report:
Gases in the atmosphere can contribute to climate change both directly and indirectly. Direct effects occur when the gas itself absorbs radiation. Indirect radiative forcing occurs when chemical transformations of the substance produce other greenhouse gases, when a gas influences the atmospheric lifetimes of other gases, and/or when a gas affects atmospheric processes that alter the radiative balance of the earth (e.g., affect cloud formation or albedo).
The Global Warming Potential (GWP) compares the ability of a greenhouse gas to trap heat in the atmosphere relative to another gas.
The GWP of a greenhouse gas is defined as the ratio of the accumulated radiative forcing within a specific time horizon caused by emitting 1 kilogram of the gas, relative to that of the reference gas CO2 (IPCC 2013). Therefore GWP-weighted emissions are provided in million metric tons of CO2 equivalent (MMT CO2 Eq.)
However, because CO2 is so much more abundant and stays in the atmosphere so much longer than other greenhouse gases, it has been the largest contributor to global warming. See here for more details.
It is also important to keep in mind that there is a lag between greenhouse gas emissions and temperature changes that we experience because much of Earth’s thermal energy (and CO2) gets stored in the ocean.
Due to a process called thermal inertia, the heat stored in the ocean will eventually be transfered to the surface of the Earth long after the gases were emitted that resulted in the increased ocean temperature.
See here for more explanation.
Furthermore, rising CO2 levels in the ocean also influence ocean acidity:
As CO2 levels rise in the ocean, the pH becomes more acidic, which makes it difficult for organisms to maintain their shells or skeletons that are made of calcium carbonate, thus making it more difficult for these organisms to survive and impacting their role in the ecosystem and food chain.
Furthermore, greenhouse gas emissions are believed to influence weather patterns as shown in this report.
Indeed, events with high levels of precipitation which can induce flooding and property damage are generally increasing around the country:
An important limitation regarding this data analysis to keep in mind is the datasets only include countries and years in which countries were reporting such information to the agencies that collected the data. Thus, the data are incomplete. For example, while we have a fairly good sense of CO2 emissions globally for later years, additional emissions were also produced by countries that are not included in the data.
In this case study we will be using data related to CO2 emissions, as well as other data that may influence, be influenced or relate to CO2 emissions. Most of our data is from Gapminder that was originally obtained from the World Bank.
In addition, we will use some data that is specific to the United States from the National Oceanic and Atmospheric Administration (NOAA), which is an agency that collects weather and climate data.
| Data | Time span | Source | Original Source | Description | Citation |
|---|---|---|---|---|---|
| CO2 emissions | 1751-2014 | Gapminder | Carbon Dioxide Information Analysis Center (CDIAC) | CO2 emissions in tonnes or metric tons (equivalent to approximately 2,204.6 pounds) per person by country | NA |
| GDP per capita (percent yearly growth) | 1801-2019 | Gapminder | World Bank | Growth Domestic Product (which is an overall measure of the health of nation’s economy) per person by country | NA |
| Energy use per person | 1960-2015 | Gapminder | World Bank | Use of primary energy before transformation to other end-use fuels, by country | NA |
| US Natural Disasters | 1980-2019 | The National Oceanic and Atmospheric Administration (NOAA) | The National Oceanic and Atmospheric Administration (NOAA) | US data about: – Droughts – Floods – Freezes – Severe Storms – Tropical Cyclones – Wildfires – Winter Storms |
NOAA National Centers for Environmental Information (NCEI) U.S. Billion-Dollar Weather and Climate Disasters (2020). https://www.ncdc.noaa.gov/billions/, DOI: 10.25921/stkw-7w73 |
| Temperature | 1895-2019 | The National Oceanic and Atmospheric Administration (NOAA) | The National Oceanic and Atmospheric Administration (NOAA) | US National yearly average temperature (in Fahrenheit) from 1895 to 2019 | NOAA National Centers for Environmental information, Climate at a Glance: National Time Series, published June 2020, retrieved on June 26, 2020 from https://www.ncdc.noaa.gov/cag/ |
To obtain the temperature data, the annual average temperatures were selected as shown in this image: 
Importantly, notice that the data we would like to use span different time periods:
| Data | Time span |
|---|---|
| CO2 emissions | 1751 to 2014 |
| GDP per capita (yearly growth) | 1801 to 2019 |
| Energy use per person | 1960 to 2015 |
| US Natural Disasters | 1980 to 2019 |
| Temperature | 1895 to 2019 |
We will explore more about this a bit later.
Question Opportunity
What concerns might arise about reliability and variation of measurement practices over time?
In our case, we downloaded the data for the files from the various sources as indicated in the table above and put them within a “raw” subdirectory of a “data” directory for our project. If you use an RStudio project, then you can use the here() function of the here package to make the path for importing this data simpler. The here package automatically starts looking for files based on where you have a .Rproj file which is created when you start a new RStudio project. We can specify that we want to look for the “yearly_co2_emissions_1000_tonnes.xlsx” file within the “raw” directory within the “data” directory within a directory where our .Rproj file is located by separating the names of these directories using commas and listing “data” first.
Click here to see more about creating new projects in RStudio.
You can create a project by going to the File menu of RStudio like so:
You can also do so by clicking the project button:
See here to learn more about using RStudio projects and here to learn more about the here package.
To read in the files that were downloaded from the various sources as indicated in the table above, we will use the read_xlsx() and read_xls() functions of the readxl package to import the data from the .xlsx and .xls files, respectively. We will also use the here() function of the here package to more easily specify the path to our files relative to the directory where the .Rproj file is located.
CO2_emissions <- readxl::read_xlsx(here("data","raw", "yearly_co2_emissions_1000_tonnes.xlsx"))
gdp_growth <- readxl::read_xlsx(here("data", "raw", "gdp_per_capita_yearly_growth.xlsx"))
energy_use <- readxl::read_xlsx(here("data", "raw", "energy_use_per_person.xlsx"))If you had trouble downloading these files, you can do so at our GitHub repo or more directly by clicking here, here, and here.
We will use the read_csv() function of the readr package to import the data from the .csv files.
However, for these files there are some lines that we would like to not import because the number of columns differ for some rows. If we don’t account for this, then we may end up importing fewer columns of the data that we would like.
In the first 5 rows shown below in the data/disasters.csv file, you can see that the first two rows does not have the same number of columns as the subsequent rows and are just (sub)titles.
To do this, we can skip rows using the skip = 2 argument of the read_csv() function.
If you had trouble downloading this file, you can do so at our GitHub repo or more directly by clicking here.
Now looking at the data/temperature.csv file, we see that the first four lines do not have the same number of columns as the subsequent lines.
We will skip importing all 4 lines by using skip = 4. We can also replace all instances of "-99" with NA using the na = "-99" argument of the read_csv() function. The “-99” needs to be in quotation marks because this argument expects characters.
Click here for an explanation about data types in R and about character strings.
There are several classes of data in R programming, meaning that certain objects will be treated or interpreted differently. Character is one of these classes. A character string is an individual data value made up of characters. This can be a paragraph, like the legend for the table, or it can be a single letter or number like the letter “a” or the number “3”. If data are of class character, than the numeric values will not be processed like a numeric value in a mathematical sense. If you want your numeric values to be interpreted that way, they need to be converted to a numeric class. The options typically used are integer (which has no decimal place) and double precision (which has a decimal place).
A variable that is a factor has a set of particular values called levels (this can be numbers or characters). Even if these are numeric, they will be interpreted as levels (i.e., as if they were characters) not as mathematical numbers. The values of a factor are assumed to have a particular ordering; by default the order is alphabetical, but this is not always the correct/intuitive ordering. You can modify the order of these levels with the forcats package.
If you had trouble downloading this file, you can do so at our GitHub repo or more directly by clicking here.
Great! now we have imported all of the data that we will need.
To allow users to skip import we will save the data as an RDA file:
If you have been following along but stopped, we could load our imported data like so:
If you skipped the data import section click here.
An RDA file (stands for R data) of the data can be found here or slightly more directly here. Download this file and then place it in your current working directory within a subdirectory called “imported” within a directory called “data” to copy and paste our code. We used an RStudio project and the here package to navigate to the file more easily.
Click here to see more about creating new projects in RStudio.
You can create a project by going to the File menu of RStudio like so:
You can also do so by clicking the project button:
See here to learn more about using RStudio projects and here to learn more about the here package.
Next, we take a look at our data that we just imported. We will need to do some data wrangling to allow us to evaluate how CO2 emissions have changed over time and how emissions may relate to energy use, GDP, etc. Let’s explore how to do that with useful functions and packages from the tidyverse.
First, let’s take a look at the CO2 data (CO2_emissions). We can use the slice_head() function of the dplyr package to see just the first rows of our data. We can specify how many rows we would like to see by using the n = argument.
We will use the %>% pipe from the magrittr package (although it is also imported by other tidyverse packages, like dplyr), which can be used to define the input for later sequential steps. This will make more sense when we have multiple sequential steps using the same data object.
# A tibble: 3 x 265
country `1751` `1752` `1753` `1754` `1755` `1756` `1757` `1758` `1759` `1760`
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Afghani… NA NA NA NA NA NA NA NA NA NA
2 Albania NA NA NA NA NA NA NA NA NA NA
3 Algeria NA NA NA NA NA NA NA NA NA NA
# … with 254 more variables: 1761 <dbl>, 1762 <dbl>, 1763 <dbl>, 1764 <dbl>,
# 1765 <dbl>, 1766 <dbl>, 1767 <dbl>, 1768 <dbl>, 1769 <dbl>, 1770 <dbl>,
# 1771 <dbl>, 1772 <dbl>, 1773 <dbl>, 1774 <dbl>, 1775 <dbl>, 1776 <dbl>,
# 1777 <dbl>, 1778 <dbl>, 1779 <dbl>, 1780 <dbl>, 1781 <dbl>, 1782 <dbl>,
# 1783 <dbl>, 1784 <dbl>, 1785 <dbl>, 1786 <dbl>, 1787 <dbl>, 1788 <dbl>,
# 1789 <dbl>, 1790 <dbl>, 1791 <dbl>, 1792 <dbl>, 1793 <dbl>, 1794 <dbl>,
# 1795 <dbl>, 1796 <dbl>, 1797 <dbl>, 1798 <dbl>, 1799 <dbl>, 1800 <dbl>,
# 1801 <dbl>, 1802 <dbl>, 1803 <dbl>, 1804 <dbl>, 1805 <dbl>, 1806 <dbl>,
# 1807 <dbl>, 1808 <dbl>, 1809 <dbl>, 1810 <dbl>, 1811 <dbl>, 1812 <dbl>,
# 1813 <dbl>, 1814 <dbl>, 1815 <dbl>, 1816 <dbl>, 1817 <dbl>, 1818 <dbl>,
# 1819 <dbl>, 1820 <dbl>, 1821 <dbl>, 1822 <dbl>, 1823 <dbl>, 1824 <dbl>,
# 1825 <dbl>, 1826 <dbl>, 1827 <dbl>, 1828 <dbl>, 1829 <dbl>, 1830 <dbl>,
# 1831 <dbl>, 1832 <dbl>, 1833 <dbl>, 1834 <dbl>, 1835 <dbl>, 1836 <dbl>,
# 1837 <dbl>, 1838 <dbl>, 1839 <dbl>, 1840 <dbl>, 1841 <dbl>, 1842 <dbl>,
# 1843 <dbl>, 1844 <dbl>, 1845 <dbl>, 1846 <dbl>, 1847 <dbl>, 1848 <dbl>,
# 1849 <dbl>, 1850 <dbl>, 1851 <dbl>, 1852 <dbl>, 1853 <dbl>, 1854 <dbl>,
# 1855 <dbl>, 1856 <dbl>, 1857 <dbl>, 1858 <dbl>, 1859 <dbl>, 1860 <dbl>, …Another useful function is slice_sample() to look at a selection of random rows using pseudorandom numbers for the index of rows to show. To continue to get the same random values or for others to get the same values, we need to set a seed first. We can do this with the set.seed() base function. We just specify a number with this function and that will allow us to get the same subset of values from the slice_sample() function. If two different people ran this code (without set.seed()), they would each see a different subset of rows. For data exploration, this isn’t a huge deal, but if we’d like separate analysts running the same code to see the same output, we will use set.seed(). If we changed set.seed(123) to set.seed(333), we would obtain a different random sample of rows.
# A tibble: 3 x 265
country `1751` `1752` `1753` `1754` `1755` `1756` `1757` `1758` `1759` `1760`
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Sri Lan… NA NA NA NA NA NA NA NA NA NA
2 Tuvalu NA NA NA NA NA NA NA NA NA NA
3 Banglad… NA NA NA NA NA NA NA NA NA NA
# … with 254 more variables: 1761 <dbl>, 1762 <dbl>, 1763 <dbl>, 1764 <dbl>,
# 1765 <dbl>, 1766 <dbl>, 1767 <dbl>, 1768 <dbl>, 1769 <dbl>, 1770 <dbl>,
# 1771 <dbl>, 1772 <dbl>, 1773 <dbl>, 1774 <dbl>, 1775 <dbl>, 1776 <dbl>,
# 1777 <dbl>, 1778 <dbl>, 1779 <dbl>, 1780 <dbl>, 1781 <dbl>, 1782 <dbl>,
# 1783 <dbl>, 1784 <dbl>, 1785 <dbl>, 1786 <dbl>, 1787 <dbl>, 1788 <dbl>,
# 1789 <dbl>, 1790 <dbl>, 1791 <dbl>, 1792 <dbl>, 1793 <dbl>, 1794 <dbl>,
# 1795 <dbl>, 1796 <dbl>, 1797 <dbl>, 1798 <dbl>, 1799 <dbl>, 1800 <dbl>,
# 1801 <dbl>, 1802 <dbl>, 1803 <dbl>, 1804 <dbl>, 1805 <dbl>, 1806 <dbl>,
# 1807 <dbl>, 1808 <dbl>, 1809 <dbl>, 1810 <dbl>, 1811 <dbl>, 1812 <dbl>,
# 1813 <dbl>, 1814 <dbl>, 1815 <dbl>, 1816 <dbl>, 1817 <dbl>, 1818 <dbl>,
# 1819 <dbl>, 1820 <dbl>, 1821 <dbl>, 1822 <dbl>, 1823 <dbl>, 1824 <dbl>,
# 1825 <dbl>, 1826 <dbl>, 1827 <dbl>, 1828 <dbl>, 1829 <dbl>, 1830 <dbl>,
# 1831 <dbl>, 1832 <dbl>, 1833 <dbl>, 1834 <dbl>, 1835 <dbl>, 1836 <dbl>,
# 1837 <dbl>, 1838 <dbl>, 1839 <dbl>, 1840 <dbl>, 1841 <dbl>, 1842 <dbl>,
# 1843 <dbl>, 1844 <dbl>, 1845 <dbl>, 1846 <dbl>, 1847 <dbl>, 1848 <dbl>,
# 1849 <dbl>, 1850 <dbl>, 1851 <dbl>, 1852 <dbl>, 1853 <dbl>, 1854 <dbl>,
# 1855 <dbl>, 1856 <dbl>, 1857 <dbl>, 1858 <dbl>, 1859 <dbl>, 1860 <dbl>, …Question Opportunity
Try setting a different seed to see the difference in the output.
OK, we see each country is represented along one row and each column contains yearly CO2 emissions. We also see that there are a lot of NA values.
We can also use the glimpse() function of the dplyr package to view our data. This allows us to see all of our variables at once. We will see a tiny bit of each variable/column with the data displayed on the right.
Rows: 192
Columns: 265
$ country <chr> "Afghanistan", "Albania", "Algeria", "Andorra", "Angola", "Ant…
$ `1751` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1752` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1753` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1754` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1755` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1756` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1757` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1758` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1759` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1760` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1761` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1762` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1763` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1764` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1765` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1766` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1767` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1768` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1769` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1770` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1771` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1772` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1773` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1774` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1775` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1776` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1777` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1778` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1779` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1780` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1781` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1782` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1783` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1784` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1785` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1786` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1787` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1788` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1789` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1790` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1791` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1792` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1793` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1794` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1795` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1796` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1797` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1798` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1799` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1800` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1801` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1802` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1803` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1804` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1805` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1806` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1807` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 169, NA, NA, NA, NA, NA, N…
$ `1808` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1809` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1810` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1811` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1812` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1813` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1814` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1815` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1816` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1817` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1818` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `1819` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 253, NA, NA, NA, NA, NA, N…
$ `1820` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 334, NA, NA, NA, NA, NA, N…
$ `1821` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 359, NA, NA, NA, NA, NA, N…
$ `1822` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 367, NA, NA, NA, NA, NA, N…
$ `1823` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 348, NA, NA, NA, NA, NA, N…
$ `1824` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 400, NA, NA, NA, NA, NA, N…
$ `1825` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 403, NA, NA, NA, NA, NA, N…
$ `1826` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 458, NA, NA, NA, NA, NA, N…
$ `1827` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 477, NA, NA, NA, NA, NA, N…
$ `1828` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 458, NA, NA, NA, NA, NA, N…
$ `1829` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 477, NA, NA, NA, NA, NA, N…
$ `1830` <dbl> NA, NA, NA, NA, NA, NA, NA, 0.032, NA, 495.000, 0.308, NA, NA,…
$ `1831` <dbl> NA, NA, NA, NA, NA, NA, NA, 3.84e-02, NA, 4.80e+02, 3.70e-01, …
$ `1832` <dbl> NA, NA, NA, NA, NA, NA, NA, 2.56e-02, NA, 5.13e+02, 2.47e-01, …
$ `1833` <dbl> NA, NA, NA, NA, NA, NA, NA, 0.032, NA, 429.000, 0.308, NA, NA,…
$ `1834` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 587, NA, NA, NA, NA, NA, N…
$ `1835` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 634, NA, NA, NA, NA, NA, N…
$ `1836` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 675, NA, NA, NA, NA, NA, N…
$ `1837` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 708, NA, NA, NA, NA, NA, N…
$ `1838` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 851, NA, NA, NA, NA, NA, N…
$ `1839` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 1060, NA, NA, NA, NA, NA, …
$ `1840` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 1170, NA, NA, NA, NA, NA, …
$ `1841` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 1320, NA, NA, NA, NA, NA, …
$ `1842` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 1460, NA, NA, NA, NA, NA, …
$ `1843` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 1270, NA, NA, NA, NA, NA, …
$ `1844` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 1600, NA, NA, NA, NA, NA, …
$ `1845` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 1800, NA, NA, NA, NA, NA, …
$ `1846` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 2120, NA, NA, NA, NA, NA, …
$ `1847` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 2080, NA, NA, NA, NA, NA, …
$ `1848` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 2340, NA, NA, NA, NA, NA, …
$ `1849` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 2260, NA, NA, NA, NA, NA, …
$ `1850` <dbl> NA, NA, NA, NA, NA, NA, NA, 0.198, NA, 2330.000, 1.910, NA, NA…
$ `1851` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 2340, NA, NA, NA, NA, NA, …
$ `1852` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 2810, NA, NA, NA, NA, NA, …
$ `1853` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 3230, NA, NA, NA, NA, NA, …
$ `1854` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 3180, NA, NA, NA, NA, NA, …
$ `1855` <dbl> NA, NA, NA, NA, NA, NA, NA, 6.01e-01, NA, 3.70e+03, 5.80e+00, …
$ `1856` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 4240, NA, NA, NA, NA, NA, …
$ `1857` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 4880, NA, NA, NA, NA, NA, …
$ `1858` <dbl> NA, NA, NA, NA, NA, NA, NA, 8.44e-01, NA, 7.25e+03, 8.14e+00, …
$ `1859` <dbl> NA, NA, NA, NA, NA, NA, NA, 8.95e-01, NA, 5.87e+03, 8.64e+00, …
$ `1860` <dbl> NA, NA, NA, NA, NA, NA, NA, 1.18, 279.00, 6150.00, 11.40, NA, …
$ `1861` <dbl> NA, NA, NA, NA, NA, NA, NA, 1.5, 510.0, 6380.0, 14.5, NA, NA, …
$ `1862` <dbl> NA, NA, NA, NA, NA, NA, NA, 1.36, 356.00, 6360.00, 13.10, NA, …
$ `1863` <dbl> NA, NA, NA, NA, NA, NA, NA, 1.42, 400.00, 5880.00, 13.70, NA, …
$ `1864` <dbl> NA, NA, NA, NA, NA, NA, NA, 1.59, 268.00, 5080.00, 15.40, NA, …
$ `1865` <dbl> NA, NA, NA, NA, NA, NA, NA, 1.52, 422.00, 5360.00, 14.70, NA, …
$ `1866` <dbl> NA, NA, NA, NA, NA, NA, NA, 4.81, 697.00, 3600.00, 46.40, NA, …
$ `1867` <dbl> NA, NA, NA, NA, NA, NA, NA, 5.52, 895.00, 4920.00, 53.20, NA, …
$ `1868` <dbl> NA, NA, NA, NA, NA, NA, NA, 4.59, 733.00, 6080.00, 44.30, NA, …
$ `1869` <dbl> NA, NA, NA, NA, NA, NA, NA, 6.23, 642.00, 6490.00, 60.10, NA, …
$ `1870` <dbl> NA, NA, NA, NA, NA, NA, NA, 6.76, 601.00, 7370.00, 65.20, NA, …
$ `1871` <dbl> NA, NA, NA, NA, NA, NA, NA, 9.12, 693.00, 10200.00, 88.00, NA,…
$ `1872` <dbl> NA, NA, NA, NA, NA, NA, NA, 9.36, 708.00, 10000.00, 90.40, NA,…
$ `1873` <dbl> NA, NA, NA, NA, NA, NA, NA, 8.79, 869.00, 10700.00, 84.80, NA,…
$ `1874` <dbl> NA, NA, NA, NA, NA, NA, NA, 10.7, 891.0, 9160.0, 103.0, NA, NA…
$ `1875` <dbl> NA, NA, NA, NA, NA, NA, NA, 12.3, 829.0, 7870.0, 119.0, NA, NA…
$ `1876` <dbl> NA, NA, NA, NA, NA, NA, NA, 15.2, 931.0, 8100.0, 147.0, NA, NA…
$ `1877` <dbl> NA, NA, NA, NA, NA, NA, NA, 15.6, 1070.0, 7290.0, 150.0, NA, N…
$ `1878` <dbl> NA, NA, NA, NA, NA, NA, NA, 20.3, 968.0, 7250.0, 196.0, NA, NA…
$ `1879` <dbl> NA, NA, NA, NA, NA, NA, NA, 20.9, 1460.0, 8870.0, 201.0, NA, N…
$ `1880` <dbl> NA, NA, NA, NA, NA, NA, NA, 24.5, 2210.0, 23700.0, 236.0, NA, …
$ `1881` <dbl> NA, NA, NA, NA, NA, NA, NA, 25.80, 1770.00, 10300.00, 249.00, …
$ `1882` <dbl> NA, NA, NA, NA, NA, NA, NA, 27.20, 2010.00, 10600.00, 262.00, …
$ `1883` <dbl> NA, NA, NA, NA, NA, NA, NA, 30.90, 2430.00, 11800.00, 298.00, …
$ `1884` <dbl> NA, NA, NA, NA, NA, NA, NA, 31.4, 2570.0, 11500.0, 303.0, NA, …
$ `1885` <dbl> NA, NA, NA, NA, NA, NA, NA, 34.20, 2910.00, 12100.00, 330.00, …
$ `1886` <dbl> NA, NA, NA, NA, NA, NA, NA, 35.10, 2890.00, 11400.00, 338.00, …
$ `1887` <dbl> NA, NA, NA, NA, NA, NA, 1090.0, 37.1, 3040.0, 12300.0, 358.0, …
$ `1888` <dbl> NA, NA, NA, NA, NA, NA, 891.0, 38.7, 3530.0, 12000.0, 373.0, N…
$ `1889` <dbl> NA, NA, NA, NA, NA, NA, 1760.0, 41.8, 3430.0, 12900.0, 403.0, …
$ `1890` <dbl> NA, NA, NA, NA, NA, NA, 1370.0, 47.3, 3550.0, 13000.0, 457.0, …
$ `1891` <dbl> NA, NA, NA, NA, NA, NA, 939.0, 52.1, 4010.0, 15000.0, 503.0, N…
$ `1892` <dbl> NA, NA, NA, NA, NA, NA, 1390.0, 55.1, 4150.0, 14500.0, 532.0, …
$ `1893` <dbl> NA, NA, NA, NA, NA, NA, 1550.0, 64.6, 3970.0, 17700.0, 624.0, …
$ `1894` <dbl> NA, NA, NA, NA, NA, NA, 1990.0, 65.8, 4360.0, 18100.0, 635.0, …
$ `1895` <dbl> NA, NA, NA, NA, NA, NA, 2270.0, 75.6, 4590.0, 20400.0, 730.0, …
$ `1896` <dbl> NA, NA, NA, NA, NA, NA, 2310, 77, 4510, 21300, 743, NA, NA, NA…
$ `1897` <dbl> NA, NA, NA, NA, NA, NA, 2080, 89, 4980, 23000, 859, NA, NA, NA…
$ `1898` <dbl> NA, NA, NA, NA, NA, NA, 2350.0, 99.9, 5620.0, 24500.0, 964.0, …
$ `1899` <dbl> NA, NA, NA, NA, NA, NA, 2920, 116, 5790, 24800, 1120, NA, NA, …
$ `1900` <dbl> NA, NA, NA, NA, NA, NA, 2070, 131, 10200, 27700, 1270, NA, NA,…
$ `1901` <dbl> NA, NA, NA, NA, NA, NA, 2490, 135, 11400, 28400, 1300, NA, NA,…
$ `1902` <dbl> NA, NA, NA, NA, NA, NA, 2820, 130, 11400, 25700, 1260, NA, NA,…
$ `1903` <dbl> NA, NA, NA, NA, NA, NA, 2860, 127, 11200, 25600, 1230, NA, NA,…
$ `1904` <dbl> NA, NA, NA, NA, NA, NA, 3800, 142, 11600, 26900, 1370, NA, NA,…
$ `1905` <dbl> NA, NA, NA, NA, NA, NA, 3990, 126, 12100, 28100, 1220, NA, NA,…
$ `1906` <dbl> NA, NA, NA, NA, NA, NA, 6260, 144, 14400, 33600, 1390, NA, NA,…
$ `1907` <dbl> NA, NA, NA, NA, NA, NA, 6260, 161, 15500, 42200, 1560, NA, NA,…
$ `1908` <dbl> NA, NA, NA, NA, NA, NA, 7620, 162, 16800, 59000, 1570, NA, NA,…
$ `1909` <dbl> NA, NA, NA, NA, NA, NA, 5940, 172, 14600, 42200, 1660, NA, NA,…
$ `1910` <dbl> NA, NA, NA, NA, NA, NA, 8910, 168, 17500, 57600, 1620, NA, NA,…
$ `1911` <dbl> NA, NA, NA, NA, NA, NA, 9950, 174, 19300, 48100, 1680, NA, NA,…
$ `1912` <dbl> NA, NA, NA, NA, NA, NA, 9490, 198, 20800, 50000, 1910, NA, NA,…
$ `1913` <dbl> NA, NA, NA, NA, NA, NA, 10200, 215, 22400, 59700, 2070, NA, NA…
$ `1914` <dbl> NA, NA, NA, NA, NA, NA, 8680, 194, 24500, 48900, 1870, NA, NA,…
$ `1915` <dbl> NA, NA, NA, NA, NA, NA, 6950, 178, 21800, 34900, 1720, NA, NA,…
$ `1916` <dbl> NA, NA, 3.67, NA, NA, NA, 4990.00, 189.00, 19300.00, 8040.00, …
$ `1917` <dbl> NA, NA, 7.33, NA, NA, NA, 2230.00, 174.00, 20800.00, 3450.00, …
$ `1918` <dbl> NA, NA, 18.3, NA, NA, NA, 2520.0, 69.5, 23000.0, 3340.0, 671.0…
$ `1919` <dbl> NA, NA, 18.3, NA, NA, NA, 3730.0, 59.4, 21800.0, 3020.0, 573.0…
$ `1920` <dbl> NA, NA, 22.0, NA, NA, NA, 5900.0, 54.2, 25800.0, 14500.0, 523.…
$ `1921` <dbl> NA, NA, 25.7, NA, NA, NA, 5540.0, 58.7, 23200.0, 19400.0, 567.…
$ `1922` <dbl> NA, NA, 25.7, NA, NA, NA, 7300.0, 71.6, 24400.0, 18600.0, 692.…
$ `1923` <dbl> NA, NA, 14.7, NA, NA, NA, 8450.0, 79.1, 24900.0, 17800.0, 764.…
$ `1924` <dbl> NA, NA, 29.3, NA, NA, NA, 11000.0, 94.3, 27100.0, 20100.0, 910…
$ `1925` <dbl> NA, NA, 33.0, NA, NA, NA, 11200.0, 93.1, 28300.0, 19000.0, 898…
$ `1926` <dbl> NA, NA, 40.3, NA, NA, NA, 11300.0, 135.0, 27900.0, 18600.0, 13…
$ `1927` <dbl> NA, NA, 58.7, NA, NA, NA, 13400.0, 168.0, 28900.0, 20100.0, 16…
$ `1928` <dbl> NA, NA, 73.30, NA, NA, NA, 12800.00, 186.00, 26300.00, 21200.0…
$ `1929` <dbl> NA, NA, 80.70, NA, NA, NA, 13100.00, 201.00, 23700.00, 24200.0…
$ `1930` <dbl> NA, NA, 84.30, NA, NA, NA, 12800.00, 273.00, 22000.00, 18900.0…
$ `1931` <dbl> NA, NA, 99.00, NA, NA, NA, 12900.00, 328.00, 19600.00, 18100.0…
$ `1932` <dbl> NA, NA, 114.00, NA, NA, NA, 13100.00, 369.00, 20400.00, 15200.…
$ `1933` <dbl> NA, 7.33, 121.00, NA, NA, NA, 13200.00, 412.00, 21600.00, 1420…
$ `1934` <dbl> NA, 7.33, 139.00, NA, NA, NA, 14300.00, 499.00, 22700.00, 1380…
$ `1935` <dbl> NA, 18.3, 132.0, NA, NA, NA, 14000.0, 565.0, 25300.0, 13900.0,…
$ `1936` <dbl> NA, 128.0, 51.3, NA, NA, NA, 15100.0, 648.0, 27100.0, 13600.0,…
$ `1937` <dbl> NA, 297.0, 69.7, NA, NA, NA, 16700.0, 662.0, 28900.0, 15300.0,…
$ `1938` <dbl> NA, 348, 33, NA, NA, NA, 16400, 699, 28100, 5790, 6750, NA, 34…
$ `1939` <dbl> NA, 433.00, 161.00, NA, NA, NA, 17400.00, 707.00, 32200.00, 63…
$ `1940` <dbl> NA, 693, 238, NA, NA, NA, 15900, 848, 29100, 7350, 8190, NA, 2…
$ `1941` <dbl> NA, 627, 312, NA, NA, NA, 14000, 745, 34600, 7980, 7190, NA, 2…
$ `1942` <dbl> NA, 744, 499, NA, NA, NA, 13500, 513, 36500, 8560, 4950, NA, 2…
$ `1943` <dbl> NA, 462, 469, NA, NA, NA, 14100, 655, 35000, 9620, 6320, NA, 2…
$ `1944` <dbl> NA, 154, 499, NA, NA, NA, 14000, 613, 34200, 9400, 5920, NA, 2…
$ `1945` <dbl> NA, 121, 616, NA, NA, NA, 13700, 649, 32700, 4570, 6270, NA, 3…
$ `1946` <dbl> NA, 484, 763, NA, NA, NA, 13700, 730, 35500, 12800, 7040, NA, …
$ `1947` <dbl> NA, 928.00, 744.00, NA, NA, NA, 14500.00, 878.00, 38000.00, 17…
$ `1948` <dbl> NA, 704.00, 803.00, NA, NA, NA, 17400.00, 935.00, 38500.00, 24…
$ `1949` <dbl> 14.70, 1020.00, 909.00, NA, NA, NA, 15400.00, 1060.00, 37700.0…
$ `1950` <dbl> 84.3, 297.0, 3790.0, NA, 187.0, NA, 30000.0, 1180.0, 54800.0, …
$ `1951` <dbl> 91.7, 403.0, 4140.0, NA, 249.0, NA, 35000.0, 1280.0, 59100.0, …
$ `1952` <dbl> 91.7, 374.0, 3890.0, NA, 312.0, NA, 36100.0, 1370.0, 60300.0, …
$ `1953` <dbl> 106.0, 414.0, 4000.0, NA, 275.0, NA, 35200.0, 1450.0, 59500.0,…
$ `1954` <dbl> 106.0, 502.0, 4160.0, NA, 348.0, NA, 36800.0, 1590.0, 67900.0,…
$ `1955` <dbl> 154.0, 664.0, 4610.0, NA, 414.0, NA, 39600.0, 1800.0, 70700.0,…
$ `1956` <dbl> 183.0, 840.0, 5000.0, NA, 502.0, NA, 44300.0, 1970.0, 73100.0,…
$ `1957` <dbl> 293.0, 1510.0, 5540.0, NA, 620.0, 22.0, 47700.0, 2160.0, 74600…
$ `1958` <dbl> 330.0, 1200.0, 5220.0, NA, 594.0, 29.3, 44200.0, 2310.0, 77700…
$ `1959` <dbl> 385.0, 1440.0, 5670.0, NA, 620.0, 29.3, 49000.0, 2430.0, 83800…
$ `1960` <dbl> 414.0, 2020.0, 6160.0, NA, 550.0, 36.7, 48800.0, 2530.0, 88200…
$ `1961` <dbl> 491.0, 2280.0, 6070.0, NA, 455.0, 47.7, 51200.0, 2600.0, 90600…
$ `1962` <dbl> 689.0, 2460.0, 5670.0, NA, 1180.0, 103.0, 53700.0, 2730.0, 949…
$ `1963` <dbl> 708.0, 2080.0, 5430.0, NA, 1150.0, 84.3, 50100.0, 2930.0, 1010…
$ `1964` <dbl> 840.0, 2020.0, 5650.0, NA, 1220.0, 91.7, 55700.0, 3120.0, 1090…
$ `1965` <dbl> 1010.0, 2170.0, 6600.0, NA, 1190.0, 150.0, 58900.0, 3310.0, 12…
$ `1966` <dbl> 1090.0, 2550.0, 8430.0, NA, 1550.0, 348.0, 63100.0, 3490.0, 12…
$ `1967` <dbl> 1280, 2680, 8440, NA, 994, 565, 65500, 3650, 129000, 40000, 35…
$ `1968` <dbl> 1220, 3070, 9060, NA, 1670, 990, 69100, 3750, 135000, 42400, 3…
$ `1969` <dbl> 942, 3250, 11300, NA, 2790, 1260, 77300, 3910, 142000, 44700, …
$ `1970` <dbl> 1.67e+03, 3.74e+03, 1.51e+04, NA, 3.58e+03, 4.62e+02, 8.27e+04…
$ `1971` <dbl> 1.90e+03, 4.35e+03, 1.87e+04, NA, 3.41e+03, 4.25e+02, 8.89e+04…
$ `1972` <dbl> 1.53e+03, 5.64e+03, 2.83e+04, NA, 4.51e+03, 3.74e+02, 9.02e+04…
$ `1973` <dbl> 1.64e+03, 5.29e+03, 3.83e+04, NA, 4.88e+03, 3.30e+02, 9.41e+04…
$ `1974` <dbl> 1.92e+03, 4.35e+03, 3.19e+04, NA, 4.87e+03, 4.29e+02, 9.56e+04…
$ `1975` <dbl> 2.13e+03, 4.59e+03, 3.20e+04, NA, 4.42e+03, 7.08e+02, 9.49e+04…
$ `1976` <dbl> 1.99e+03, 4.95e+03, 3.92e+04, NA, 3.29e+03, 4.03e+02, 9.98e+04…
$ `1977` <dbl> 2.39e+03, 5.72e+03, 4.19e+04, NA, 3.53e+03, 4.66e+02, 1.01e+05…
$ `1978` <dbl> 2160, 6490, 62500, NA, 5410, 491, 103000, 5810, 202000, 57500,…
$ `1979` <dbl> 2240, 7590, 45600, NA, 5500, 407, 111000, 5850, 205000, 61600,…
$ `1980` <dbl> 1760, 5170, 66500, NA, 5350, 143, 109000, 6080, 221000, 52300,…
$ `1981` <dbl> 1980.0, 7340.0, 46400.0, NA, 5280.0, 106.0, 102000.0, 5970.0, …
$ `1982` <dbl> 2100, 7310, 39300, NA, 4650, 293, 103000, 6080, 234000, 53900,…
$ `1983` <dbl> 2520.0, 7630.0, 52600.0, NA, 5120.0, 84.3, 105000.0, 6170.0, 2…
$ `1984` <dbl> 2830.0, 7830.0, 71100.0, NA, 5010.0, 147.0, 107000.0, 6230.0, …
$ `1985` <dbl> 3510.0, 7880.0, 72800.0, NA, 4700.0, 249.0, 101000.0, 6710.0, …
$ `1986` <dbl> 3140, 8060, 76300, NA, 4660, 249, 104000, 6730, 240000, 54100,…
$ `1987` <dbl> 3120, 7440, 84100, NA, 5820, 275, 115000, 7020, 256000, 57700,…
$ `1988` <dbl> 2870, 7330, 83900, NA, 5130, 286, 121000, 7210, 261000, 53300,…
$ `1989` <dbl> 2780.0, 8980.0, 80000.0, NA, 5010.0, 286.0, 117000.0, 7060.0, …
$ `1990` <dbl> 2610, 5520, 77000, 407, 5120, 282, 112000, 6620, 264000, 57700…
$ `1991` <dbl> 2440, 4290, 79000, 407, 5090, 268, 117000, 6380, 261000, 61600…
$ `1992` <dbl> 1390, 2520, 80100, 407, 5200, 264, 121000, 5830, 268000, 56700…
$ `1993` <dbl> 1350, 2340, 82200, 411, 5780, 271, 118000, 2560, 277000, 57100…
$ `1994` <dbl> 1290, 1930, 86400, 407, 3890, 268, 122000, 2710, 278000, 57100…
$ `1995` <dbl> 1240, 2090, 95300, 425, 11000, 275, 128000, 3410, 282000, 5980…
$ `1996` <dbl> 1180, 2020, 97100, 455, 10500, 293, 135000, 2560, 302000, 6320…
$ `1997` <dbl> 1100, 1540, 87300, 466, 7380, 308, 138000, 3230, 306000, 62700…
$ `1998` <dbl> 1040, 1750, 107000, 491, 7310, 319, 140000, 3360, 317000, 6370…
$ `1999` <dbl> 821, 2980, 92000, 513, 9160, 330, 147000, 3010, 325000, 61900,…
$ `2000` <dbl> 774, 3020, 87900, 524, 9540, 345, 142000, 3470, 329000, 62300,…
$ `2001` <dbl> 818, 3220, 84200, 524, 9730, 348, 134000, 3540, 325000, 65900,…
$ `2002` <dbl> 1070, 3750, 89900, 532, 12700, 370, 125000, 3040, 341000, 6710…
$ `2003` <dbl> 1200, 4290, 91600, 535, 9060, 403, 135000, 3430, 336000, 72200…
$ `2004` <dbl> 950, 4170, 88500, 561, 18800, 422, 158000, 3640, 343000, 72400…
$ `2005` <dbl> 1330, 4250, 107000, 576, 19200, 429, 162000, 4350, 350000, 742…
$ `2006` <dbl> 1650, 3900, 101000, 546, 22300, 444, 175000, 4380, 365000, 722…
$ `2007` <dbl> 2270, 3930, 109000, 539, 25200, 469, 175000, 5060, 372000, 697…
$ `2008` <dbl> 4210, 4370, 110000, 539, 25700, 480, 189000, 5560, 386000, 690…
$ `2009` <dbl> 6770, 4380, 121000, 517, 27800, 510, 180000, 4360, 395000, 627…
$ `2010` <dbl> 8460, 4600, 119000, 517, 29100, 524, 188000, 4220, 391000, 675…
$ `2011` <dbl> 12200, 5240, 121000, 491, 30300, 513, 192000, 4920, 392000, 65…
$ `2012` <dbl> 10800, 4910, 130000, 488, 33400, 524, 192000, 5690, 388000, 62…
$ `2013` <dbl> 10000, 5060, 134000, 477, 32600, 524, 190000, 5500, 372000, 62…
$ `2014` <dbl> 9810, 5720, 145000, 462, 34800, 532, 204000, 5530, 361000, 587…We can also see that we have a large tibble.
[1] "tbl_df" "tbl" "data.frame"This is the object that is created when we read in the data with readr. A tibble (or tbl_df) is the tidyverse version of a data.frame object. Similar to data.frame, it is a table with variable information arranged as columns, and individual observations arranged as rows. However some nice differences are they do not change variable names or data types and they give more messages when something is wrong (e.g. when a variable does not exist), which forces the analyst to confront problems earlier. Tibbles also give us information about the class of each variable.
For example the country variable is made up of character (abbreviated as chr) values.
# A tibble: 192 x 1
country
<chr>
1 Afghanistan
2 Albania
3 Algeria
4 Andorra
5 Angola
6 Antigua and Barbuda
7 Argentina
8 Armenia
9 Australia
10 Austria
# … with 182 more rowsWe see that we have 192 rows different country variables and CO2 emission values for 264 different years (from 1751 to 2014).
[1] "country" "1751" "1752" "1753" "1754" "1755" "1756"
[8] "1757" "1758" "1759" "1760" "1761" "1762" "1763"
[15] "1764" "1765" "1766" "1767" "1768" "1769" "1770"
[22] "1771" "1772" "1773" "1774" "1775" "1776" "1777"
[29] "1778" "1779" "1780" "1781" "1782" "1783" "1784"
[36] "1785" "1786" "1787" "1788" "1789" "1790" "1791"
[43] "1792" "1793" "1794" "1795" "1796" "1797" "1798"
[50] "1799" "1800" "1801" "1802" "1803" "1804" "1805"
[57] "1806" "1807" "1808" "1809" "1810" "1811" "1812"
[64] "1813" "1814" "1815" "1816" "1817" "1818" "1819"
[71] "1820" "1821" "1822" "1823" "1824" "1825" "1826"
[78] "1827" "1828" "1829" "1830" "1831" "1832" "1833"
[85] "1834" "1835" "1836" "1837" "1838" "1839" "1840"
[92] "1841" "1842" "1843" "1844" "1845" "1846" "1847"
[99] "1848" "1849" "1850" "1851" "1852" "1853" "1854"
[106] "1855" "1856" "1857" "1858" "1859" "1860" "1861"
[113] "1862" "1863" "1864" "1865" "1866" "1867" "1868"
[120] "1869" "1870" "1871" "1872" "1873" "1874" "1875"
[127] "1876" "1877" "1878" "1879" "1880" "1881" "1882"
[134] "1883" "1884" "1885" "1886" "1887" "1888" "1889"
[141] "1890" "1891" "1892" "1893" "1894" "1895" "1896"
[148] "1897" "1898" "1899" "1900" "1901" "1902" "1903"
[155] "1904" "1905" "1906" "1907" "1908" "1909" "1910"
[162] "1911" "1912" "1913" "1914" "1915" "1916" "1917"
[169] "1918" "1919" "1920" "1921" "1922" "1923" "1924"
[176] "1925" "1926" "1927" "1928" "1929" "1930" "1931"
[183] "1932" "1933" "1934" "1935" "1936" "1937" "1938"
[190] "1939" "1940" "1941" "1942" "1943" "1944" "1945"
[197] "1946" "1947" "1948" "1949" "1950" "1951" "1952"
[204] "1953" "1954" "1955" "1956" "1957" "1958" "1959"
[211] "1960" "1961" "1962" "1963" "1964" "1965" "1966"
[218] "1967" "1968" "1969" "1970" "1971" "1972" "1973"
[225] "1974" "1975" "1976" "1977" "1978" "1979" "1980"
[232] "1981" "1982" "1983" "1984" "1985" "1986" "1987"
[239] "1988" "1989" "1990" "1991" "1992" "1993" "1994"
[246] "1995" "1996" "1997" "1998" "1999" "2000" "2001"
[253] "2002" "2003" "2004" "2005" "2006" "2007" "2008"
[260] "2009" "2010" "2011" "2012" "2013" "2014" Recall, the values are emissions in metric tons, also called tonnes. Scrolling through the glimpse() function above, we can also see that there are fewer NA values for later years.
In this next code chunk, we will introduce the %<>% operator from the magrittr package. This allows us to use our CO2_emissions data and reassign it to a modified version at the same time. Let’s modify CO2_emissions to make it more usable for making visualizations. Specifically, we will use the pivot_longer() function of the dplyr package to convert our data into what is called “long” format. This is also sometimes referred to as “narrow” format.
This means that we will have more rows and fewer columns than our current format.
Right now our data is in what is called “wide” format. In wide format, each variable is listed as its own column. In contrast, in long format, variables maybe collapsed into a column that identifies the variables and a column of values. See here for more information about the difference between the two formats.
We want to collapse all of the values for the emission data across the different individual year variables into one new Emissions variable. We will identify what year they are from by creating a new Year variable. The cols = argument allows us to specify which columns we want to pivot (or not pivot) to create these new columns. We want to keep our country data as an ID variable so we will exclude it using the - sign, by default all other columns will be used.
CO2_emissions %<>%
pivot_longer(cols = -country,
names_to = "Year",
values_to = "Emissions")
set.seed(123)
CO2_emissions %>%
slice_sample(n = 6)# A tibble: 6 x 3
country Year Emissions
<chr> <chr> <dbl>
1 Bahamas 1832 NA
2 Montenegro 1843 NA
3 Mongolia 1892 NA
4 Samoa 1791 NA
5 Azerbaijan 1867 53.2
6 Timor-Leste 2010 235 Question Opportunity
Think a moment about what the dimensions of the CO2_emissions tibble are now and why? How would you check this?
Hint : Checking has something to do with a unique aspect about tibbles.
Let’s say we also want to rename the country variable to be capitalized. To do this, we can use the rename() function of the dplyr package to rename this variable. When renaming variables the syntax is new-name = old-name, where the new name is listed first before the =.
You may also note that the Year variable is currently of class type character. We would like to change it to be numeric. To do this we will use the mutate() function, which is also part of the dplyr package. This function allows us to create and modify variables. We will also use this function to create a variable called Label which will have "CO2 Emissions (Metric Tons)" as the value for every row, to be used when we create plots later.
CO2_emissions %<>%
dplyr::rename(Country = country) %>%
dplyr::mutate(Year = as.numeric(Year),
Label = "CO2 Emissions (Metric Tons)")Now let’s take a look to see how our data has changed:
# A tibble: 6 x 4
Country Year Emissions Label
<chr> <dbl> <dbl> <chr>
1 Bahamas 1832 NA CO2 Emissions (Metric Tons)
2 Montenegro 1843 NA CO2 Emissions (Metric Tons)
3 Mongolia 1892 NA CO2 Emissions (Metric Tons)
4 Samoa 1791 NA CO2 Emissions (Metric Tons)
5 Azerbaijan 1867 53.2 CO2 Emissions (Metric Tons)
6 Timor-Leste 2010 235 CO2 Emissions (Metric Tons)Great, we can see that now the Year variable is of class double (abbreviated dbl), which is a numeric class.
Now, let’s take a look at the Country variable to check if there is anything unexpected. We will use the distinct() function of the dplyr package to view the unique values only. Finally, we use the pull() function of the dplyr package to extract the values from the column (this is similar to using the $ base R syntax e.g. CO2_emission$Country).
[1] "Afghanistan" "Albania"
[3] "Algeria" "Andorra"
[5] "Angola" "Antigua and Barbuda"
[7] "Argentina" "Armenia"
[9] "Australia" "Austria"
[11] "Azerbaijan" "Bahamas"
[13] "Bahrain" "Bangladesh"
[15] "Barbados" "Belarus"
[17] "Belgium" "Belize"
[19] "Benin" "Bhutan"
[21] "Bolivia" "Bosnia and Herzegovina"
[23] "Botswana" "Brazil"
[25] "Brunei" "Bulgaria"
[27] "Burkina Faso" "Burundi"
[29] "Cambodia" "Cameroon"
[31] "Canada" "Cape Verde"
[33] "Central African Republic" "Chad"
[35] "Chile" "China"
[37] "Colombia" "Comoros"
[39] "Congo, Dem. Rep." "Congo, Rep."
[41] "Costa Rica" "Cote d'Ivoire"
[43] "Croatia" "Cuba"
[45] "Cyprus" "Czech Republic"
[47] "Denmark" "Djibouti"
[49] "Dominica" "Dominican Republic"
[51] "Ecuador" "Egypt"
[53] "El Salvador" "Equatorial Guinea"
[55] "Eritrea" "Estonia"
[57] "Ethiopia" "Fiji"
[59] "Finland" "France"
[61] "Gabon" "Gambia"
[63] "Georgia" "Germany"
[65] "Ghana" "Greece"
[67] "Grenada" "Guatemala"
[69] "Guinea" "Guinea-Bissau"
[71] "Guyana" "Haiti"
[73] "Honduras" "Hungary"
[75] "Iceland" "India"
[77] "Indonesia" "Iran"
[79] "Iraq" "Ireland"
[81] "Israel" "Italy"
[83] "Jamaica" "Japan"
[85] "Jordan" "Kazakhstan"
[87] "Kenya" "Kiribati"
[89] "Kuwait" "Kyrgyz Republic"
[91] "Lao" "Latvia"
[93] "Lebanon" "Lesotho"
[95] "Liberia" "Libya"
[97] "Liechtenstein" "Lithuania"
[99] "Luxembourg" "Macedonia, FYR"
[101] "Madagascar" "Malawi"
[103] "Malaysia" "Maldives"
[105] "Mali" "Malta"
[107] "Marshall Islands" "Mauritania"
[109] "Mauritius" "Mexico"
[111] "Micronesia, Fed. Sts." "Moldova"
[113] "Mongolia" "Montenegro"
[115] "Morocco" "Mozambique"
[117] "Myanmar" "Namibia"
[119] "Nauru" "Nepal"
[121] "Netherlands" "New Zealand"
[123] "Nicaragua" "Niger"
[125] "Nigeria" "North Korea"
[127] "Norway" "Oman"
[129] "Pakistan" "Palau"
[131] "Palestine" "Panama"
[133] "Papua New Guinea" "Paraguay"
[135] "Peru" "Philippines"
[137] "Poland" "Portugal"
[139] "Qatar" "Romania"
[141] "Russia" "Rwanda"
[143] "Samoa" "Sao Tome and Principe"
[145] "Saudi Arabia" "Senegal"
[147] "Serbia" "Seychelles"
[149] "Sierra Leone" "Singapore"
[151] "Slovak Republic" "Slovenia"
[153] "Solomon Islands" "Somalia"
[155] "South Africa" "South Korea"
[157] "South Sudan" "Spain"
[159] "Sri Lanka" "St. Kitts and Nevis"
[161] "St. Lucia" "St. Vincent and the Grenadines"
[163] "Sudan" "Suriname"
[165] "Swaziland" "Sweden"
[167] "Switzerland" "Syria"
[169] "Tajikistan" "Tanzania"
[171] "Thailand" "Timor-Leste"
[173] "Togo" "Tonga"
[175] "Trinidad and Tobago" "Tunisia"
[177] "Turkey" "Turkmenistan"
[179] "Tuvalu" "Uganda"
[181] "Ukraine" "United Arab Emirates"
[183] "United Kingdom" "United States"
[185] "Uruguay" "Uzbekistan"
[187] "Vanuatu" "Venezuela"
[189] "Vietnam" "Yemen"
[191] "Zambia" "Zimbabwe" These all look as expected!
Let’s take a look at the next dataset (gdp_growth) that we imported.
# A tibble: 3 x 220
country `1801` `1802` `1803` `1804` `1805` `1806` `1807`
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Afghanistan NA NA NA NA NA NA NA
2 Albania 0.104 0.104 0.104 0.104 0.104 0.104 0.104
3 Algeria -0.00247 -0.00247 -0.00247 -0.00247 -0.00247 -0.00247 -0.00247
# … with 212 more variables: 1808 <dbl>, 1809 <dbl>, 1810 <dbl>, 1811 <dbl>,
# 1812 <dbl>, 1813 <dbl>, 1814 <dbl>, 1815 <dbl>, 1816 <dbl>, 1817 <dbl>,
# 1818 <dbl>, 1819 <dbl>, 1820 <dbl>, 1821 <dbl>, 1822 <dbl>, 1823 <dbl>,
# 1824 <dbl>, 1825 <dbl>, 1826 <dbl>, 1827 <dbl>, 1828 <dbl>, 1829 <dbl>,
# 1830 <dbl>, 1831 <dbl>, 1832 <dbl>, 1833 <dbl>, 1834 <dbl>, 1835 <dbl>,
# 1836 <dbl>, 1837 <dbl>, 1838 <dbl>, 1839 <dbl>, 1840 <dbl>, 1841 <dbl>,
# 1842 <dbl>, 1843 <dbl>, 1844 <dbl>, 1845 <dbl>, 1846 <dbl>, 1847 <dbl>,
# 1848 <dbl>, 1849 <dbl>, 1850 <dbl>, 1851 <dbl>, 1852 <dbl>, 1853 <dbl>,
# 1854 <dbl>, 1855 <dbl>, 1856 <dbl>, 1857 <dbl>, 1858 <dbl>, 1859 <dbl>,
# 1860 <dbl>, 1861 <dbl>, 1862 <dbl>, 1863 <dbl>, 1864 <dbl>, 1865 <dbl>,
# 1866 <dbl>, 1867 <dbl>, 1868 <dbl>, 1869 <dbl>, 1870 <dbl>, 1871 <dbl>,
# 1872 <dbl>, 1873 <dbl>, 1874 <dbl>, 1875 <dbl>, 1876 <dbl>, 1877 <dbl>,
# 1878 <dbl>, 1879 <dbl>, 1880 <dbl>, 1881 <dbl>, 1882 <dbl>, 1883 <dbl>,
# 1884 <dbl>, 1885 <dbl>, 1886 <dbl>, 1887 <dbl>, 1888 <dbl>, 1889 <dbl>,
# 1890 <dbl>, 1891 <dbl>, 1892 <dbl>, 1893 <dbl>, 1894 <dbl>, 1895 <dbl>,
# 1896 <dbl>, 1897 <dbl>, 1898 <dbl>, 1899 <dbl>, 1900 <dbl>, 1901 <dbl>,
# 1902 <dbl>, 1903 <dbl>, 1904 <dbl>, 1905 <dbl>, 1906 <dbl>, 1907 <dbl>, …How many rows and columns are there are there? We can easily check by using the base dim() function, which evaluates the dimensions of an object.
[1] 194 220Interesting, it’s 194 rows (as opposed to 50688 above). We will deal with this and other differences in the sets of countries a bit later on. There are also 220 columns with a country column and a set of columns corresponding to different years.
[1] "country" "1801" "1802" "1803" "1804" "1805" "1806"
[8] "1807" "1808" "1809" "1810" "1811" "1812" "1813"
[15] "1814" "1815" "1816" "1817" "1818" "1819" "1820"
[22] "1821" "1822" "1823" "1824" "1825" "1826" "1827"
[29] "1828" "1829" "1830" "1831" "1832" "1833" "1834"
[36] "1835" "1836" "1837" "1838" "1839" "1840" "1841"
[43] "1842" "1843" "1844" "1845" "1846" "1847" "1848"
[50] "1849" "1850" "1851" "1852" "1853" "1854" "1855"
[57] "1856" "1857" "1858" "1859" "1860" "1861" "1862"
[64] "1863" "1864" "1865" "1866" "1867" "1868" "1869"
[71] "1870" "1871" "1872" "1873" "1874" "1875" "1876"
[78] "1877" "1878" "1879" "1880" "1881" "1882" "1883"
[85] "1884" "1885" "1886" "1887" "1888" "1889" "1890"
[92] "1891" "1892" "1893" "1894" "1895" "1896" "1897"
[99] "1898" "1899" "1900" "1901" "1902" "1903" "1904"
[106] "1905" "1906" "1907" "1908" "1909" "1910" "1911"
[113] "1912" "1913" "1914" "1915" "1916" "1917" "1918"
[120] "1919" "1920" "1921" "1922" "1923" "1924" "1925"
[127] "1926" "1927" "1928" "1929" "1930" "1931" "1932"
[134] "1933" "1934" "1935" "1936" "1937" "1938" "1939"
[141] "1940" "1941" "1942" "1943" "1944" "1945" "1946"
[148] "1947" "1948" "1949" "1950" "1951" "1952" "1953"
[155] "1954" "1955" "1956" "1957" "1958" "1959" "1960"
[162] "1961" "1962" "1963" "1964" "1965" "1966" "1967"
[169] "1968" "1969" "1970" "1971" "1972" "1973" "1974"
[176] "1975" "1976" "1977" "1978" "1979" "1980" "1981"
[183] "1982" "1983" "1984" "1985" "1986" "1987" "1988"
[190] "1989" "1990" "1991" "1992" "1993" "1994" "1995"
[197] "1996" "1997" "1998" "1999" "2000" "2001" "2002"
[204] "2003" "2004" "2005" "2006" "2007" "2008" "2009"
[211] "2010" "2011" "2012" "2013" "2014" "2015" "2016"
[218] "2017" "2018" "2019" Yes, no other columns in this dataset.
Next, we will use the pivot_longer() to transform the data to long format, similar to what we did in the previous section.
We will also again change the country variable to be Country by using the rename() function, and we will make the Year variable numeric using the mutate() function.
Question Opportunity
Using what you just learned about pivot_longer(), rename(), and mutate() and without scrolling up, try to come up with the code to do the wrangling for this data.
Click here to reveal the code.
Now let’s see how this data has changed:
# A tibble: 6 x 4
Country Year GDP Label
<chr> <dbl> <dbl> <chr>
1 Afghanistan 1801 NA GDP Growth/Capita (%)
2 Afghanistan 1802 NA GDP Growth/Capita (%)
3 Afghanistan 1803 NA GDP Growth/Capita (%)
4 Afghanistan 1804 NA GDP Growth/Capita (%)
5 Afghanistan 1805 NA GDP Growth/Capita (%)
6 Afghanistan 1806 NA GDP Growth/Capita (%)# A tibble: 219 x 2
Year n
* <dbl> <int>
1 1801 194
2 1802 194
3 1803 194
4 1804 194
5 1805 194
6 1806 194
7 1807 194
8 1808 194
9 1809 194
10 1810 194
# … with 209 more rowsAgain let’s check that the Country variable only contains values we would expect.
[1] "Afghanistan" "Albania"
[3] "Algeria" "Andorra"
[5] "Angola" "Antigua and Barbuda"
[7] "Argentina" "Armenia"
[9] "Australia" "Austria"
[11] "Azerbaijan" "Bahamas"
[13] "Bahrain" "Bangladesh"
[15] "Barbados" "Belarus"
[17] "Belgium" "Belize"
[19] "Benin" "Bhutan"
[21] "Bolivia" "Bosnia and Herzegovina"
[23] "Botswana" "Brazil"
[25] "Brunei" "Bulgaria"
[27] "Burkina Faso" "Burundi"
[29] "Cambodia" "Cameroon"
[31] "Canada" "Cape Verde"
[33] "Central African Republic" "Chad"
[35] "Chile" "China"
[37] "Colombia" "Comoros"
[39] "Congo, Dem. Rep." "Congo, Rep."
[41] "Costa Rica" "Cote d'Ivoire"
[43] "Croatia" "Cuba"
[45] "Cyprus" "Czech Republic"
[47] "Denmark" "Djibouti"
[49] "Dominica" "Dominican Republic"
[51] "Ecuador" "Egypt"
[53] "El Salvador" "Equatorial Guinea"
[55] "Eritrea" "Estonia"
[57] "Ethiopia" "Fiji"
[59] "Finland" "France"
[61] "Gabon" "Gambia"
[63] "Georgia" "Germany"
[65] "Ghana" "Greece"
[67] "Grenada" "Guatemala"
[69] "Guinea" "Guinea-Bissau"
[71] "Guyana" "Haiti"
[73] "Honduras" "Hungary"
[75] "Iceland" "India"
[77] "Indonesia" "Iran"
[79] "Iraq" "Ireland"
[81] "Israel" "Italy"
[83] "Jamaica" "Japan"
[85] "Jordan" "Kazakhstan"
[87] "Kenya" "Kiribati"
[89] "Kuwait" "Kyrgyz Republic"
[91] "Lao" "Latvia"
[93] "Lebanon" "Lesotho"
[95] "Liberia" "Libya"
[97] "Liechtenstein" "Lithuania"
[99] "Luxembourg" "Macedonia, FYR"
[101] "Madagascar" "Malawi"
[103] "Malaysia" "Maldives"
[105] "Mali" "Malta"
[107] "Marshall Islands" "Mauritania"
[109] "Mauritius" "Mexico"
[111] "Micronesia, Fed. Sts." "Moldova"
[113] "Monaco" "Mongolia"
[115] "Montenegro" "Morocco"
[117] "Mozambique" "Myanmar"
[119] "Namibia" "Nauru"
[121] "Nepal" "Netherlands"
[123] "New Zealand" "Nicaragua"
[125] "Niger" "Nigeria"
[127] "North Korea" "Norway"
[129] "Oman" "Pakistan"
[131] "Palau" "Palestine"
[133] "Panama" "Papua New Guinea"
[135] "Paraguay" "Peru"
[137] "Philippines" "Poland"
[139] "Portugal" "Qatar"
[141] "Romania" "Russia"
[143] "Rwanda" "Samoa"
[145] "San Marino" "Sao Tome and Principe"
[147] "Saudi Arabia" "Senegal"
[149] "Serbia" "Seychelles"
[151] "Sierra Leone" "Singapore"
[153] "Slovak Republic" "Slovenia"
[155] "Solomon Islands" "Somalia"
[157] "South Africa" "South Korea"
[159] "South Sudan" "Spain"
[161] "Sri Lanka" "St. Kitts and Nevis"
[163] "St. Lucia" "St. Vincent and the Grenadines"
[165] "Sudan" "Suriname"
[167] "Swaziland" "Sweden"
[169] "Switzerland" "Syria"
[171] "Tajikistan" "Tanzania"
[173] "Thailand" "Timor-Leste"
[175] "Togo" "Tonga"
[177] "Trinidad and Tobago" "Tunisia"
[179] "Turkey" "Turkmenistan"
[181] "Tuvalu" "Uganda"
[183] "Ukraine" "United Arab Emirates"
[185] "United Kingdom" "United States"
[187] "Uruguay" "Uzbekistan"
[189] "Vanuatu" "Venezuela"
[191] "Vietnam" "Yemen"
[193] "Zambia" "Zimbabwe" Also looks good!
Now let’s take a look at the energy use per person data (energy_use) using slice_head() and glimpse().
# A tibble: 3 x 57
country `1960` `1961` `1962` `1963` `1964` `1965` `1966` `1967` `1968` `1969`
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Albania NA NA NA NA NA NA NA NA NA NA
2 Algeria NA NA NA NA NA NA NA NA NA NA
3 Angola NA NA NA NA NA NA NA NA NA NA
# … with 46 more variables: 1970 <dbl>, 1971 <dbl>, 1972 <dbl>, 1973 <dbl>,
# 1974 <dbl>, 1975 <dbl>, 1976 <dbl>, 1977 <dbl>, 1978 <dbl>, 1979 <dbl>,
# 1980 <dbl>, 1981 <dbl>, 1982 <dbl>, 1983 <dbl>, 1984 <dbl>, 1985 <dbl>,
# 1986 <dbl>, 1987 <dbl>, 1988 <dbl>, 1989 <dbl>, 1990 <dbl>, 1991 <dbl>,
# 1992 <dbl>, 1993 <dbl>, 1994 <dbl>, 1995 <dbl>, 1996 <dbl>, 1997 <dbl>,
# 1998 <dbl>, 1999 <dbl>, 2000 <dbl>, 2001 <dbl>, 2002 <dbl>, 2003 <dbl>,
# 2004 <dbl>, 2005 <dbl>, 2006 <dbl>, 2007 <dbl>, 2008 <dbl>, 2009 <dbl>,
# 2010 <dbl>, 2011 <dbl>, 2012 <dbl>, 2013 <dbl>, 2014 <dbl>, 2015 <dbl>Rows: 169
Columns: 57
$ country <chr> "Albania", "Algeria", "Angola", "Antigua and Barbuda", "Argent…
$ `1960` <dbl> NA, NA, NA, NA, NA, NA, 3060, 1550, NA, NA, NA, NA, NA, NA, 25…
$ `1961` <dbl> NA, NA, NA, NA, NA, NA, 3120, 1550, NA, NA, NA, NA, NA, NA, 25…
$ `1962` <dbl> NA, NA, NA, NA, NA, NA, 3170, 1680, NA, NA, NA, NA, NA, NA, 28…
$ `1963` <dbl> NA, NA, NA, NA, NA, NA, 3280, 1820, NA, NA, NA, NA, NA, NA, 30…
$ `1964` <dbl> NA, NA, NA, NA, NA, NA, 3350, 1860, NA, NA, NA, NA, NA, NA, 30…
$ `1965` <dbl> NA, NA, NA, NA, NA, NA, 3460, 1850, NA, NA, NA, NA, NA, NA, 31…
$ `1966` <dbl> NA, NA, NA, NA, NA, NA, 3550, 1900, NA, NA, NA, NA, NA, NA, 30…
$ `1967` <dbl> NA, NA, NA, NA, NA, NA, 3690, 1920, NA, NA, NA, NA, NA, NA, 31…
$ `1968` <dbl> NA, NA, NA, NA, NA, NA, 3760, 2050, NA, NA, NA, NA, NA, NA, 35…
$ `1969` <dbl> NA, NA, NA, NA, NA, NA, 3790, 2180, NA, NA, NA, NA, NA, NA, 38…
$ `1970` <dbl> NA, NA, NA, NA, NA, NA, 4060, 2420, NA, NA, NA, NA, NA, NA, 41…
$ `1971` <dbl> 785.0, 232.0, 556.0, NA, 1380.0, NA, 3990.0, 2510.0, NA, NA, 6…
$ `1972` <dbl> 866.0, 261.0, 584.0, NA, 1380.0, NA, 4040.0, 2630.0, NA, NA, 5…
$ `1973` <dbl> 763.0, 305.0, 568.0, NA, 1410.0, NA, 4260.0, 2830.0, NA, NA, 8…
$ `1974` <dbl> 777.0, 319.0, 565.0, NA, 1420.0, NA, 4290.0, 2730.0, NA, NA, 9…
$ `1975` <dbl> 827.0, 330.0, 536.0, NA, 1380.0, NA, 4350.0, 2650.0, NA, NA, 8…
$ `1976` <dbl> 891, 367, 515, NA, 1400, NA, 4410, 2870, NA, NA, 9580, 98, NA,…
$ `1977` <dbl> 924.0, 399.0, 494.0, NA, 1420.0, NA, 4670.0, 2800.0, NA, NA, 8…
$ `1978` <dbl> 1010.0, 477.0, 527.0, NA, 1430.0, NA, 4630.0, 2890.0, NA, NA, …
$ `1979` <dbl> 864.0, 586.0, 518.0, NA, 1480.0, NA, 4680.0, 3140.0, NA, NA, 8…
$ `1980` <dbl> 1150, 579, 511, NA, 1490, NA, 4740, 3070, NA, NA, 7790, 103, N…
$ `1981` <dbl> 989, 611, 497, NA, 1430, NA, 4690, 2900, NA, NA, 8300, 102, NA…
$ `1982` <dbl> 967, 771, 473, NA, 1420, NA, 4820, 2830, NA, NA, 9070, 105, NA…
$ `1983` <dbl> 1000, 808, 469, NA, 1420, NA, 4560, 2840, NA, NA, 8500, 105, N…
$ `1984` <dbl> 1020, 776, 458, NA, 1450, NA, 4650, 2950, NA, NA, 8830, 104, N…
$ `1985` <dbl> 917, 786, 470, NA, 1360, NA, 4600, 3050, NA, NA, 9920, 107, NA…
$ `1986` <dbl> 964, 862, 462, NA, 1420, NA, 4620, 3060, NA, NA, 10300, 111, N…
$ `1987` <dbl> 922, 828, 461, NA, 1480, NA, 4770, 3170, NA, NA, 9520, 107, NA…
$ `1988` <dbl> 928, 850, 467, NA, 1500, NA, 4700, 3200, NA, NA, 10500, 114, N…
$ `1989` <dbl> 896, 820, 465, NA, 1440, NA, 5000, 3140, NA, NA, 10200, 117, N…
$ `1990` <dbl> 813, 856, 483, 1480, 1410, 2180, 5060, 3240, 3170, 2520, 10600…
$ `1991` <dbl> 573, 884, 480, NA, 1430, 2320, 4930, 3420, 3090, NA, 10100, 11…
$ `1992` <dbl> 418, 884, 467, NA, 1480, 1200, 4960, 3250, 2460, NA, 10800, 11…
$ `1993` <dbl> 412, 868, 468, NA, 1470, 652, 5150, 3260, 2180, NA, 11100, 123…
$ `1994` <dbl> 441, 819, 459, NA, 1540, 420, 5090, 3230, 1950, NA, 11600, 126…
$ `1995` <dbl> 417, 839, 445, NA, 1540, 511, 5130, 3370, 1810, NA, 11400, 134…
$ `1996` <dbl> 448, 798, 445, NA, 1580, 562, 5390, 3580, 1510, NA, 11100, 132…
$ `1997` <dbl> 385, 805, 443, NA, 1610, 594, 5470, 3550, 1440, NA, 12200, 135…
$ `1998` <dbl> 427, 821, 430, NA, 1650, 610, 5550, 3610, 1490, NA, 12400, 138…
$ `1999` <dbl> 576, 864, 439, NA, 1660, 594, 5610, 3590, 1370, NA, 11900, 137…
$ `2000` <dbl> 580, 866, 437, NA, 1660, 656, 5640, 3570, 1400, NA, 12000, 139…
$ `2001` <dbl> 597, 856, 442, NA, 1560, 657, 5450, 3760, 1410, NA, 11700, 149…
$ `2002` <dbl> 660, 904, 447, NA, 1500, 618, 5570, 3770, 1410, NA, 11500, 150…
$ `2003` <dbl> 648, 949, 466, NA, 1590, 657, 5570, 3970, 1480, NA, 11600, 155…
$ `2004` <dbl> 715, 948, 462, 1530, 1720, 698, 5600, 4010, 1540, 2060, 10900,…
$ `2005` <dbl> 720, 974, 431, 1530, 1710, 843, 5560, 4090, 1600, 2110, 11700,…
$ `2006` <dbl> 707, 1030, 456, 1580, 1840, 865, 5710, 4080, 1560, 2100, 11600…
$ `2007` <dbl> 680, 1070, 470, 1600, 1850, 973, 5870, 4020, 1410, 2070, 11200…
$ `2008` <dbl> 711, 1070, 491, NA, 1920, 1030, 5960, 4030, 1520, NA, 11300, 1…
$ `2009` <dbl> 732, 1150, 514, NA, 1850, 904, 5860, 3800, 1330, NA, 10300, 18…
$ `2010` <dbl> 729, 1110, 521, NA, 1910, 863, 5790, 4050, 1280, NA, 10200, 20…
$ `2011` <dbl> 765, 1140, 522, NA, 1930, 944, 5750, 3920, 1370, NA, 9910, 206…
$ `2012` <dbl> 688, 1220, 553, NA, 1920, 1030, 5570, 3890, 1470, NA, 9660, 21…
$ `2013` <dbl> 801, 1240, 534, NA, 1950, 1000, 5460, 3920, 1470, NA, 10400, 2…
$ `2014` <dbl> 808, 1320, 545, NA, 2020, 1020, 5330, 3760, 1500, NA, 10600, 2…
$ `2015` <dbl> NA, NA, NA, NA, NA, NA, 5480, 3800, NA, NA, NA, NA, NA, NA, 46…Looks like we have 169 rows and 57 columns where we have a country column and again a set of years. To wrangle the energy_use data, we will again convert the data to long format, rename some variables, and mutate the Year data to be numeric.
Question Opportunity
Again try to come up with the code on your own to wrangle the data.
Click here to reveal the code.
# A tibble: 3 x 4
Country Year Energy Label
<chr> <dbl> <dbl> <chr>
1 Dominican Republic 2014 734 Energy Use (kg, oil-eq./capita)
2 Ecuador 2006 671 Energy Use (kg, oil-eq./capita)
3 Turkey 1997 1170 Energy Use (kg, oil-eq./capita)Now we will check the Country variable:
[1] "Albania" "Algeria"
[3] "Angola" "Antigua and Barbuda"
[5] "Argentina" "Armenia"
[7] "Australia" "Austria"
[9] "Azerbaijan" "Bahamas"
[11] "Bahrain" "Bangladesh"
[13] "Barbados" "Belarus"
[15] "Belgium" "Belize"
[17] "Benin" "Bhutan"
[19] "Bolivia" "Bosnia and Herzegovina"
[21] "Botswana" "Brazil"
[23] "Brunei" "Bulgaria"
[25] "Cambodia" "Cameroon"
[27] "Canada" "Cape Verde"
[29] "Chile" "China"
[31] "Colombia" "Comoros"
[33] "Congo, Dem. Rep." "Congo, Rep."
[35] "Costa Rica" "Cote d'Ivoire"
[37] "Croatia" "Cuba"
[39] "Cyprus" "Czech Republic"
[41] "Denmark" "Djibouti"
[43] "Dominica" "Dominican Republic"
[45] "Ecuador" "Egypt"
[47] "El Salvador" "Equatorial Guinea"
[49] "Eritrea" "Estonia"
[51] "Ethiopia" "Fiji"
[53] "Finland" "France"
[55] "Gabon" "Gambia"
[57] "Georgia" "Germany"
[59] "Ghana" "Greece"
[61] "Grenada" "Guatemala"
[63] "Guinea-Bissau" "Guyana"
[65] "Haiti" "Honduras"
[67] "Hungary" "Iceland"
[69] "India" "Indonesia"
[71] "Iran" "Iraq"
[73] "Ireland" "Israel"
[75] "Italy" "Jamaica"
[77] "Japan" "Jordan"
[79] "Kazakhstan" "Kenya"
[81] "Kiribati" "Kuwait"
[83] "Kyrgyz Republic" "Latvia"
[85] "Lebanon" "Lesotho"
[87] "Libya" "Lithuania"
[89] "Luxembourg" "Macedonia, FYR"
[91] "Malaysia" "Maldives"
[93] "Malta" "Marshall Islands"
[95] "Mauritius" "Mexico"
[97] "Moldova" "Mongolia"
[99] "Montenegro" "Morocco"
[101] "Mozambique" "Myanmar"
[103] "Namibia" "Nepal"
[105] "Netherlands" "New Zealand"
[107] "Nicaragua" "Niger"
[109] "Nigeria" "North Korea"
[111] "Norway" "Oman"
[113] "Pakistan" "Palau"
[115] "Panama" "Paraguay"
[117] "Peru" "Philippines"
[119] "Poland" "Portugal"
[121] "Qatar" "Romania"
[123] "Russia" "Samoa"
[125] "Sao Tome and Principe" "Saudi Arabia"
[127] "Senegal" "Serbia"
[129] "Seychelles" "Singapore"
[131] "Slovak Republic" "Slovenia"
[133] "Solomon Islands" "South Africa"
[135] "South Korea" "South Sudan"
[137] "Spain" "Sri Lanka"
[139] "St. Kitts and Nevis" "St. Lucia"
[141] "St. Vincent and the Grenadines" "Sudan"
[143] "Suriname" "Swaziland"
[145] "Sweden" "Switzerland"
[147] "Syria" "Tajikistan"
[149] "Tanzania" "Thailand"
[151] "Timor-Leste" "Togo"
[153] "Tonga" "Trinidad and Tobago"
[155] "Tunisia" "Turkey"
[157] "Turkmenistan" "Ukraine"
[159] "United Arab Emirates" "United Kingdom"
[161] "United States" "Uruguay"
[163] "Uzbekistan" "Vanuatu"
[165] "Venezuela" "Vietnam"
[167] "Yemen" "Zambia"
[169] "Zimbabwe" Looks good!
Now we will take a look at the US data about disasters and temperature.
First, we consider the disasters that have occurred in the US.
# A tibble: 40 x 57
Year `Drought Count` `Drought Cost` `Drought Lower 75` `Drought Upper 75`
<dbl> <dbl> <dbl> <dbl> <dbl>
1 1980 1 33.2 26.4 39.6
2 1981 0 0 0 0
3 1982 0 0 0 0
4 1983 1 7.8 5.5 9
5 1984 0 0 0 0
6 1985 0 0 0 0
7 1986 1 4.2 3.5 5
8 1987 0 0 0 0
9 1988 1 44.4 33.8 54
10 1989 1 6.4 5.6 7.4
# … with 30 more rows, and 52 more variables: Drought Lower 90 <dbl>,
# Drought Upper 90 <dbl>, Drought Lower 95 <dbl>, Drought Upper 95 <dbl>,
# Flooding Count <dbl>, Flooding Cost <dbl>, Flooding Lower 75 <dbl>,
# Flooding Upper 75 <dbl>, Flooding Lower 90 <dbl>, Flooding Upper 90 <dbl>,
# Flooding Lower 95 <dbl>, Flooding Upper 95 <dbl>, Freeze Count <dbl>,
# Freeze Cost <dbl>, Freeze Lower 75 <dbl>, Freeze Upper 75 <dbl>,
# Freeze Lower 90 <dbl>, Freeze Upper 90 <dbl>, Freeze Lower 95 <dbl>,
# Freeze Upper 95 <dbl>, Severe Storm Count <dbl>, Severe Storm Cost <dbl>,
# Severe Storm Lower 75 <dbl>, Severe Storm Upper 75 <dbl>,
# Severe Storm Lower 90 <dbl>, Severe Storm Upper 90 <dbl>,
# Severe Storm Lower 95 <dbl>, Severe Storm Upper 95 <dbl>,
# Tropical Cyclone Count <dbl>, Tropical Cyclone Cost <dbl>,
# Tropical Cyclone Lower 75 <dbl>, Tropical Cyclone Upper 75 <dbl>,
# Tropical Cyclone Lower 90 <dbl>, Tropical Cyclone Upper 90 <dbl>,
# Tropical Cyclone Lower 95 <dbl>, Tropical Cyclone Upper 95 <dbl>,
# Wildfire Count <dbl>, Wildfire Cost <dbl>, Wildfire Lower 75 <dbl>,
# Wildfire Upper 75 <dbl>, Wildfire Lower 90 <dbl>, Wildfire Upper 90 <dbl>,
# Wildfire Lower 95 <dbl>, Wildfire Upper 95 <dbl>, Winter Storm Count <dbl>,
# Winter Storm Cost <dbl>, Winter Storm Lower 75 <dbl>,
# Winter Storm Upper 75 <dbl>, Winter Storm Lower 90 <dbl>,
# Winter Storm Upper 90 <dbl>, Winter Storm Lower 95 <dbl>,
# Winter Storm Upper 95 <dbl>We are specifically interested in the Year and the variables that contain the word "Count". The other variables represent an estimate of the economic cost in billions of dollars, as well as the upper and lower bounds for simulations used to estimate the economic cost, which show the level of uncertainty in these estimates (at three different levels of confidence) as the true cost is unknown. See here for more information about the data. For this analysis, we will focus just on the number of disasters that occurred each year.
We will select our variables of interest using the select() and contains() functions in the dplyr package. Since we are selecting for variables with the word "Count" we need to use quotation marks around it.
Selecting for the variable Year does not require quotes because it is the full name of one of the existing variables.
# A tibble: 6 x 8
Year `Drought Count` `Flooding Count` `Freeze Count` `Severe Storm Count`
<dbl> <dbl> <dbl> <dbl> <dbl>
1 1980 1 1 0 0
2 1981 0 0 1 1
3 1982 0 0 0 2
4 1983 1 2 1 0
5 1984 0 0 0 2
6 1985 0 0 1 0
# … with 3 more variables: Tropical Cyclone Count <dbl>, Wildfire Count <dbl>,
# Winter Storm Count <dbl>Now we want to create a new variable that will be the sum of all the different types of disasters for each year.
We don’t want to include the Year variable in our sum, so we can exclude it using the select function. To perform the sum for each year, we can use the base rowSums() function.
[1] 3 2 3 5 2 5 2 0 1 5 3 4 7 5 6 5 4 3 10 5 4 2 4 7 5
[26] 6 7 5 12 7 6 16 11 9 8 10 15 16 14 14We could then add this to our us_diaster tibble like so using the bind_cols function of the dplyr package:
# A tibble: 40 x 9
Year `Drought Count` `Flooding Count` `Freeze Count` `Severe Storm Count`
<dbl> <dbl> <dbl> <dbl> <dbl>
1 1980 1 1 0 0
2 1981 0 0 1 1
3 1982 0 0 0 2
4 1983 1 2 1 0
5 1984 0 0 0 2
6 1985 0 0 1 0
7 1986 1 0 0 1
8 1987 0 0 0 0
9 1988 1 0 0 0
10 1989 1 0 1 1
# … with 30 more rows, and 4 more variables: Tropical Cyclone Count <dbl>,
# Wildfire Count <dbl>, Winter Storm Count <dbl>, Disaters <dbl>However, we can actually create and add this new variable directly to the us_disaster tibble by using the mutate() function of dplyr and using the . notation.
We need to use the . notation to indicate that we are using the data that we already used as input (on the left side of the pipe) to our mutate() function (on the right side of our pipe), which in this case is the entire us_disaster tibble for our select() function. The output from the select() function will be used for the rowSums() function.
Rows: 40
Columns: 9
$ Year <dbl> 1980, 1981, 1982, 1983, 1984, 1985, 1986, 198…
$ `Drought Count` <dbl> 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, …
$ `Flooding Count` <dbl> 1, 0, 0, 2, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, …
$ `Freeze Count` <dbl> 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, …
$ `Severe Storm Count` <dbl> 0, 1, 2, 0, 2, 0, 1, 0, 0, 1, 1, 1, 4, 1, 1, …
$ `Tropical Cyclone Count` <dbl> 1, 0, 0, 1, 0, 3, 0, 0, 0, 1, 0, 1, 2, 0, 1, …
$ `Wildfire Count` <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, …
$ `Winter Storm Count` <dbl> 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 2, …
$ Disasters <dbl> 3, 2, 3, 5, 2, 5, 2, 0, 1, 5, 3, 4, 7, 5, 6, …Great, now we are going to remove some of these variables and just keep the variables of interest using select().
We are also going to add a new variable called Country to indicate that this data is from the United States. Again this will create a new variable where every value is United States.
us_disaster %<>%
dplyr::select(Year, Disasters) %>%
mutate(Country = "United States") %>%
pivot_longer(cols = c(-Country, -Year),
names_to = "Indicator",
values_to = "Value") %>%
mutate(Label = "Number of Disasters")
us_disaster %>%
slice_head(n = 6)# A tibble: 6 x 5
Year Country Indicator Value Label
<dbl> <chr> <chr> <dbl> <chr>
1 1980 United States Disasters 3 Number of Disasters
2 1981 United States Disasters 2 Number of Disasters
3 1982 United States Disasters 3 Number of Disasters
4 1983 United States Disasters 5 Number of Disasters
5 1984 United States Disasters 2 Number of Disasters
6 1985 United States Disasters 5 Number of DisastersGreat, this looks good now.
Question Opportunity
This dataset was slightly different from the other datasets and therefore required slightly different wrangling. Why was it necessary to exclude the Year variable from the pivot_longer() function? What would happen if we did not exclude Year?
Next, we consider the temperature in the US over time.
# A tibble: 6 x 3
Date Value Anomaly
<dbl> <dbl> <dbl>
1 189512 50.3 -1.68
2 189612 52.0 -0.03
3 189712 51.6 -0.46
4 189812 51.4 -0.59
5 189912 51.0 -1.01
6 190012 52.8 0.75So a few things need to be fixed here.
First, the Date column looks a bit strange. The format of the numbers look like the year followed by the number 12 (representing 12 months).
We want to change this to only keep the first 4 characters in the Date variable string values.
However, first let’s make sure that indeed all of the Date variables are 6 characters long and that they all end with the number 12.
We can use a couple of functions in the stringr package to do this. This package is used for working with character strings. The str_length() function can be used to check the length of each value, while the str_ends() function can be used to check that all the values end with "12".
Let’s start with the str_length() function. These functions in the stringr package require a character vector. Thus we need to pull the values for the Date variable first using the pull() function of the dplyr package.
[1] 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
[38] 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
[75] 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
[112] 6 6 6 6 6 6 6 6 6 6 6 6 6 6Great! It looks like all of the values are 6 characters long.
Now let’s check that they all end with "12". We just need to specify what pattern to look for.
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[106] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[121] TRUE TRUE TRUE TRUE TRUEGreat! Since all of the values are TRUE we know that all of the values in the Date variable end with "12".
It’s a good idea to always check that your data is as you expect.
Now we can use the str_sub() function of the stringr package to remove the "12" from each Date value.
We just need to indicate the start and stop characters.
In this case the start would be 1 and the 4th character would be where we want to stop, so we would use start = 1, stop = 4. We can do this inside of the mutate() function to modify the Date variable. In doing so, we will not need to use pull() to pull the values for the Date variable.
# A tibble: 125 x 3
Date Value Anomaly
<chr> <dbl> <dbl>
1 1895 50.3 -1.68
2 1896 52.0 -0.03
3 1897 51.6 -0.46
4 1898 51.4 -0.59
5 1899 51.0 -1.01
6 1900 52.8 0.75
7 1901 51.9 -0.15
8 1902 51.6 -0.43
9 1903 50.6 -1.4
10 1904 51.2 -0.86
# … with 115 more rowsWe also want to remove the Anomaly variable, which is an indicator of how different the national average temperature for that year was from the average temperature from 1901-2000 which was 52.02°F.
Then, we also want to create a Country variable. We will also change the name of the Date variable to Year so that it will be consistent with our other datasets. We also also want the Year to be numeric. We can accomplish both renaming and changing to numeric by using the mutate() function.
We also want to create an Indicator variable so that we can later tell what data the values in this tibble represent if we combine it with other tibbles and a Label variable, so that we will have informative labels if we make a plot with this data later.
Finally, we remove the Date variable and also order the columns just like the other us data using the select() function.
us_temperature %<>%
dplyr::select(-Anomaly) %>%
mutate(Country = "United States",
Year = as.numeric(Date),
Indicator = "Temperature",
Label = "Temperature (Fahrenheit)") %>%
select(Year, Country, Indicator, Value, Label)
us_temperature %>%
slice_head(n = 6)# A tibble: 6 x 5
Year Country Indicator Value Label
<dbl> <chr> <chr> <dbl> <chr>
1 1895 United States Temperature 50.3 Temperature (Fahrenheit)
2 1896 United States Temperature 52.0 Temperature (Fahrenheit)
3 1897 United States Temperature 51.6 Temperature (Fahrenheit)
4 1898 United States Temperature 51.4 Temperature (Fahrenheit)
5 1899 United States Temperature 51.0 Temperature (Fahrenheit)
6 1900 United States Temperature 52.8 Temperature (Fahrenheit)Now that we have wrangled the individual datasets, we are ready to put everything together. Specifically, we will join the individual datasets into one tibble using *_join() functions available in the dplyr package.
Before we begin though, we will need to make sure that there is at least one variable/column that has the same name across all datasets to be joined. Such variables with common names are called keys for joining your data.
These are the by="x1" arguments below where x1 is the name of the column in both the a and b datasets that we will join together.
There are several types of *_join() functions to consider. The full_join() function keeps all rows from both tibbles that are being joined and adds NA values as necessary if there are values within the key for either of the tibbles that is is not in the key of the other tibble.
We use the full_join() function as we have different time spans for each dataset and we would like to retain as much data as possible.
The full_join() function will simply create NA values for any of the years that are not in one of the datasets.
First, we check using the base summary() function that there are column names that are consistent in each dataset that we wish to combine.
Country Year Emissions Label
Length:50688 Min. :1751 Min. : 0 Length:50688
Class :character 1st Qu.:1817 1st Qu.: 550 Class :character
Mode :character Median :1882 Median : 4390 Mode :character
Mean :1882 Mean : 83808
3rd Qu.:1948 3rd Qu.: 31925
Max. :2014 Max. :10300000
NA's :33772 Country Year GDP Label
Length:42486 Min. :1801 Min. :-67.500 Length:42486
Class :character 1st Qu.:1855 1st Qu.: 0.133 Class :character
Mode :character Median :1910 Median : 0.633 Mode :character
Mean :1910 Mean : 1.302
3rd Qu.:1965 3rd Qu.: 2.160
Max. :2019 Max. :145.000
NA's :2392 Country Year Energy Label
Length:9464 Min. :1960 Min. : 9.58 Length:9464
Class :character 1st Qu.:1974 1st Qu.: 505.75 Class :character
Mode :character Median :1988 Median : 1185.00 Mode :character
Mean :1988 Mean : 2238.82
3rd Qu.:2001 3rd Qu.: 3030.00
Max. :2015 Max. :22000.00
NA's :3544 Question Opportunity
What variable or variables might we want to use to join our data by?
Click here to see an explanation for what variable or variables to join by after you have thought about it.
The Country, and Year variables are present in all of the datasets with values that overlap. Although Label is also present in the datasets, the values do not overlap. We can see that the minimum and maximum year is different for nearly all the datasets.
Next, we need to specify what columns/variables we will be joining by using the by = argument in the full_join() function, (recall that this variable is called the “key”)
data_wide <- CO2_emissions %>%
full_join(gdp_growth, by = c("Country", "Year", "Label")) %>%
full_join(energy_use, by = c("Country", "Year", "Label"))
set.seed(123)
data_wide %>%
slice_sample(n = 6)# A tibble: 6 x 6
Country Year Emissions Label GDP Energy
<chr> <dbl> <dbl> <chr> <dbl> <dbl>
1 Angola 1899 NA GDP Growth/Capita (%) 0.425 NA
2 Central African Repu… 1974 NA GDP Growth/Capita (%) -4.92 NA
3 Bahamas 1832 NA CO2 Emissions (Metric Ton… NA NA
4 Montenegro 1843 NA CO2 Emissions (Metric Ton… NA NA
5 Croatia 2015 NA Energy Use (kg, oil-eq./c… NA NA
6 Israel 1885 NA GDP Growth/Capita (%) 0.582 NA Click here to see an explanation for another option that works well for large numbers of tibbles
We can also do the same thing by using the reduce() function of the purrr package. This takes a list of elements (which can be tibbles) and then applies a function that requires two inputs iteratively using the first pair of elements and creating a single element and then applying the function again to the output element and the next listed element and so on and so forth.
For example we will use a list of tibbles and the full_join() function which requires two tibbles to combine. This will first combine CO2_emissions and gdp_growth and then take the resulting joined tibble and combine this with the energy_use tibble.
You can see that this is a great option if you have many datasets to combine!
data_wide <-
list(CO2_emissions, gdp_growth, energy_use) %>%
reduce(full_join, by = c("Country", "Year", "Label"))
set.seed(123)
data_wide %>%
slice_sample(n = 6)# A tibble: 6 x 6
Country Year Emissions Label GDP Energy
<chr> <dbl> <dbl> <chr> <dbl> <dbl>
1 Angola 1899 NA GDP Growth/Capita (%) 0.425 NA
2 Central African Repu… 1974 NA GDP Growth/Capita (%) -4.92 NA
3 Bahamas 1832 NA CO2 Emissions (Metric Ton… NA NA
4 Montenegro 1843 NA CO2 Emissions (Metric Ton… NA NA
5 Croatia 2015 NA Energy Use (kg, oil-eq./c… NA NA
6 Israel 1885 NA GDP Growth/Capita (%) 0.582 NARows: 102,638
Columns: 6
$ Country <chr> "Afghanistan", "Afghanistan", "Afghanistan", "Afghanistan", …
$ Year <dbl> 1751, 1752, 1753, 1754, 1755, 1756, 1757, 1758, 1759, 1760, …
$ Emissions <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ Label <chr> "CO2 Emissions (Metric Tons)", "CO2 Emissions (Metric Tons)"…
$ GDP <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ Energy <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …Nice, looks good!
We will also make a long version of this data, where we will create an new variable called Indicator that will indicate what dataset the data came from.
Question Opportunity
Try to come up with the code to do this.
Click here to reveal the code.
# A tibble: 6 x 5
Country Year Label Indicator Value
<chr> <dbl> <chr> <chr> <dbl>
1 Ethiopia 1829 GDP Growth/Capita (%) GDP 0.0197
2 Tanzania 1820 CO2 Emissions (Metric Tons) Energy NA
3 South Sudan 1907 CO2 Emissions (Metric Tons) GDP NA
4 Bangladesh 1931 GDP Growth/Capita (%) GDP -2
5 Mongolia 1805 GDP Growth/Capita (%) Emissions NA
6 Spain 1805 CO2 Emissions (Metric Tons) Emissions NA We will now combine this data with the US data about disasters and temperatures.
# A tibble: 6 x 5
Year Country Indicator Value Label
<dbl> <chr> <chr> <dbl> <chr>
1 1980 United States Disasters 3 Number of Disasters
2 1981 United States Disasters 2 Number of Disasters
3 1982 United States Disasters 3 Number of Disasters
4 1983 United States Disasters 5 Number of Disasters
5 1984 United States Disasters 2 Number of Disasters
6 1985 United States Disasters 5 Number of Disasters# A tibble: 6 x 5
Year Country Indicator Value Label
<dbl> <chr> <chr> <dbl> <chr>
1 1895 United States Temperature 50.3 Temperature (Fahrenheit)
2 1896 United States Temperature 52.0 Temperature (Fahrenheit)
3 1897 United States Temperature 51.6 Temperature (Fahrenheit)
4 1898 United States Temperature 51.4 Temperature (Fahrenheit)
5 1899 United States Temperature 51.0 Temperature (Fahrenheit)
6 1900 United States Temperature 52.8 Temperature (Fahrenheit)We will now use the bind_rows() function of the dplyr package which will just append the us_temperature data and the us_disaster data after the data_long data.
data_long <-
list(data_long, us_disaster, us_temperature) %>%
bind_rows() %>%
mutate(Country = as.factor(Country))We also converted the Country column to a factor in the last line of the code chunk.
We can check the top and bottom of the new data_long tibble to see that our us_temperature data is at the bottom. To see the end of our tibble we can use slice_tail() function of the dplyr package.
# A tibble: 6 x 5
Country Year Label Indicator Value
<fct> <dbl> <chr> <chr> <dbl>
1 Afghanistan 1751 CO2 Emissions (Metric Tons) Emissions NA
2 Afghanistan 1751 CO2 Emissions (Metric Tons) GDP NA
3 Afghanistan 1751 CO2 Emissions (Metric Tons) Energy NA
4 Afghanistan 1752 CO2 Emissions (Metric Tons) Emissions NA
5 Afghanistan 1752 CO2 Emissions (Metric Tons) GDP NA
6 Afghanistan 1752 CO2 Emissions (Metric Tons) Energy NA# A tibble: 6 x 5
Country Year Label Indicator Value
<fct> <dbl> <chr> <chr> <dbl>
1 United States 2014 Temperature (Fahrenheit) Temperature 52.5
2 United States 2015 Temperature (Fahrenheit) Temperature 54.4
3 United States 2016 Temperature (Fahrenheit) Temperature 54.9
4 United States 2017 Temperature (Fahrenheit) Temperature 54.6
5 United States 2018 Temperature (Fahrenheit) Temperature 53.5
6 United States 2019 Temperature (Fahrenheit) Temperature 52.7# A tibble: 10 x 5
Country Year Label Indicator Value
<fct> <dbl> <chr> <chr> <dbl>
1 Ethiopia 1829 GDP Growth/Capita (%) GDP 0.0197
2 Tanzania 1820 CO2 Emissions (Metric Tons) Energy NA
3 South Sudan 1907 CO2 Emissions (Metric Tons) GDP NA
4 Bangladesh 1931 GDP Growth/Capita (%) GDP -2
5 Mongolia 1805 GDP Growth/Capita (%) Emissions NA
6 Spain 1805 CO2 Emissions (Metric Tons) Emissions NA
7 Germany 1858 GDP Growth/Capita (%) Emissions NA
8 Fiji 1837 CO2 Emissions (Metric Tons) GDP NA
9 Jamaica 1823 CO2 Emissions (Metric Tons) Emissions NA
10 Iceland 1926 CO2 Emissions (Metric Tons) Emissions NA Click here for details about the difference between
full_join() and bind_rows()
The difference between this function and the full_join() function is that the bind_rows() function will essentially just append each dataset to each other, whereas the full_join() function collapses data that is comparable. Here, you will see an example of what the data would have been like for data_wide if we had made it using bind_rows() and if full_join() had been used but was not joined by the Label variable. Since the Label variable has unique values for each type of Indicator, this causes the full_join() result to be the same as bind_rows().
Let’s consider an example and look at the values for China in the year of 1980.
First we will use the bind_rows() function which automatically creates NA values for any variable that is missing from a data object that is added by combining the data object with another that contains that missing variable using this function.
data_wide_br <-
list(CO2_emissions, gdp_growth, energy_use) %>%
bind_rows()
data_wide_br %>%
filter(Country == "China",
Year == 1980)# A tibble: 3 x 6
Country Year Emissions Label GDP Energy
<chr> <dbl> <dbl> <chr> <dbl> <dbl>
1 China 1980 1470000 CO2 Emissions (Metric Tons) NA NA
2 China 1980 NA GDP Growth/Capita (%) 2.16 NA
3 China 1980 NA Energy Use (kg, oil-eq./capita) NA 609We see that we have three rows of data.
Now we will use the full_join() function two ways. First we will combine by Country and Year and Label.
data_wide_fj_label <-
list(CO2_emissions, gdp_growth, energy_use) %>%
reduce(full_join, by = c("Country", "Year", "Label"))
data_wide_fj_label %>%
filter(Country == "China", Year == "1980")# A tibble: 3 x 6
Country Year Emissions Label GDP Energy
<chr> <dbl> <dbl> <chr> <dbl> <dbl>
1 China 1980 1470000 CO2 Emissions (Metric Tons) NA NA
2 China 1980 NA GDP Growth/Capita (%) 2.16 NA
3 China 1980 NA Energy Use (kg, oil-eq./capita) NA 609Again we have 3 rows of data. The data produced by bind_rows() and full_join() is identical (which we can check by using the setequal() function of the dplyr package) and has the same dimensions (which we can check by using the dim() base function).
[1] 102638 6[1] 102638 6[1] TRUEHowever, now we will join by only Country and Year:
data_wide_fj <-
list(CO2_emissions, gdp_growth, energy_use) %>%
reduce(full_join, by = c("Country", "Year"))
data_wide_fj %>%
filter(Country == "China", Year == "1980")# A tibble: 1 x 8
Country Year Emissions Label.x GDP Label.y Energy Label
<chr> <dbl> <dbl> <chr> <dbl> <chr> <dbl> <chr>
1 China 1980 1470000 CO2 Emissions… 2.16 GDP Growth/… 609 Energy Use (…Now we see that we have only a single row. The data that corresponds to the same year and country has been collapsed into a single but wider row.
This is something to keep in mind when you are wrangling your data. The choice of what function to use and how should depend on how you want the data to be after you combine the different sources of data together.
We have a few more things to do before we leave the data wrangling section.
We will create a new variable called Region that will indicate if the data is about the United States or a different country based on the values in the Country variable. To do this, we will use the case_when() function of the dplyr package.
For example, if the Country variable is equal to "United States" the value for the new variable will also be "United States", where as if the Country variable is not equal to "United States" but is some other character string value, such as "Afghanistan", then the value for the new variable will be "Rest of the World". We can specify that something is not equal by using the != operator.
The new values for the new variable Region are indicated after the specific conditional statements by using the ~ symbol.
data_long %<>%
mutate(Region = case_when(Country == "United States" ~ "United States",
Country != "United States" ~ "Rest of the World"))
data_long %>%
arrange(Country) %>%
slice_head(n = 6)# A tibble: 6 x 6
Country Year Label Indicator Value Region
<fct> <dbl> <chr> <chr> <dbl> <chr>
1 Afghanistan 1751 CO2 Emissions (Metric Tons) Emissions NA Rest of the Wor…
2 Afghanistan 1751 CO2 Emissions (Metric Tons) GDP NA Rest of the Wor…
3 Afghanistan 1751 CO2 Emissions (Metric Tons) Energy NA Rest of the Wor…
4 Afghanistan 1752 CO2 Emissions (Metric Tons) Emissions NA Rest of the Wor…
5 Afghanistan 1752 CO2 Emissions (Metric Tons) GDP NA Rest of the Wor…
6 Afghanistan 1752 CO2 Emissions (Metric Tons) Energy NA Rest of the Wor…We can also remove rows for countries with NA values using the drop_na() function of the tidyr package to drop all years with missing data.
You can see that by removing the NA values the data for Afghanistan starts at 1949 instead of 1751.
# A tibble: 6 x 6
Country Year Label Indicator Value Region
<fct> <dbl> <chr> <chr> <dbl> <chr>
1 Afghanistan 1949 CO2 Emissions (Metric Tons) Emissions 14.7 Rest of the Wor…
2 Afghanistan 1950 CO2 Emissions (Metric Tons) Emissions 84.3 Rest of the Wor…
3 Afghanistan 1951 CO2 Emissions (Metric Tons) Emissions 91.7 Rest of the Wor…
4 Afghanistan 1952 CO2 Emissions (Metric Tons) Emissions 91.7 Rest of the Wor…
5 Afghanistan 1953 CO2 Emissions (Metric Tons) Emissions 106 Rest of the Wor…
6 Afghanistan 1954 CO2 Emissions (Metric Tons) Emissions 106 Rest of the Wor…Question Opportunity
Using only data from the US, calculate the first and last year that a value was reported for each variable (e.g. CO2 emissions, energy use, etc).
Hint : Use the group_by() and summarize() functions in the dplyr package.
To allow users to skip import and wrangling we will save the data as an RDA file as well as a CSV file as this is often useful to send our data to collaborators. We will save this in a “wrangled” subdirectory of our “data” directory of our working directory.
We have created a dataset called wide_format that is in “wide” format. Write some code to convert it to “long” format. (We want 3 columns with names “ID”, “Variable”, and “Value”.)
If you have been following along but stopped, we could load our imported data like so:
If you skipped the data wrangling section click here.
An RDA file (stands for R data) of the data can be found here or slightly more directly here. Download this file and then place it in your current working directory within a subdirectory called “wrangled” within a directory called “data” to copy and paste our code. We used an RStudio project and the here package to navigate to the file more easily.
Click here to see more about creating new projects in RStudio.
You can create a project by going to the File menu of RStudio like so:
You can also do so by clicking the project button:
See here to learn more about using RStudio projects and here to learn more about the here package.
Now we will create some simple plots to examine the CO2 emissions over time using the ggplot2 package.
As you may have already seen, there are many functions available in base R that can create plots (e.g. plot(), boxplot()). Others include: hist(), qqplot(), etc. These functions are great because they come with a basic installation of R and can be quite powerful when you need a quick visualization of something when you are exploring data.
We are choosing to introduce ggplot2 because, in our opinion, it is one of the simplest ways for beginners to create relatively complicated plots that are intuitive and aesthetically pleasing.
ggplot2 R packageThe reasons ggplot2 is generally intuitive for beginners is the use of grammar of graphics or the gg in ggplot2. The idea is that you can construct many sentences by learning just a few nouns, adjectives, and verbs. There are specific “words” that we will introduce, and once you are comfortable with them, you will be able to create (or “write”) hundreds of different plots.
The critical part to making graphics using ggplot2 is the data needs to be in a tidy format. Given that we have just spent time putting our data in tidy format, we are primed to take advantage of all that ggplot2 has to offer!
We will show how it is easy to pipe tidy data (output) as input to other functions that create plots. This all works because we are working within the tidyverse.
ggplot() function?As explained by Hadley Wickham:
the grammar tells us that a statistical graphic is a mapping from data to aesthetic attributes (colour, shape, size) of geometric objects (points, lines, bars). The plot may also contain statistical transformations of the data and is drawn on a specific coordinates system.
ggplot2 Terminologyx and y variable names)geom_point(), geom_bar(), geom_line(), geom_histogram()The function aes() is an aesthetic mapping function inside the ggplot() object. We use this function to specify plot attributes (e.g. x and y variable names) that will not change as we add more layers.
Anything that goes in the ggplot() object becomes a global setting. From there, we use the geom objects to add more layers to the base ggplot() object. These will define what we are interested in illustrating using the data.
Let’s start by plotting the CO2 emissions over time. Because our dataset contains other variables, we first need to filter our data to only include the CO2 emissions data by using the filter() function of the dplyr package. To use this function we need to specify what value (e.g. Emissions) we want for a given variable or column (e.g. Indicator).
In this case, we filter to keep all rows where the Indicator variable is equal to the word Emissions. Notice that this needs to be in quotes, while the variable name does not.
# A tibble: 16,916 x 6
Country Year Label Indicator Value Region
<fct> <dbl> <chr> <chr> <dbl> <chr>
1 Afghanistan 1949 CO2 Emissions (Metric Ton… Emissions 14.7 Rest of the Wor…
2 Afghanistan 1950 CO2 Emissions (Metric Ton… Emissions 84.3 Rest of the Wor…
3 Afghanistan 1951 CO2 Emissions (Metric Ton… Emissions 91.7 Rest of the Wor…
4 Afghanistan 1952 CO2 Emissions (Metric Ton… Emissions 91.7 Rest of the Wor…
5 Afghanistan 1953 CO2 Emissions (Metric Ton… Emissions 106 Rest of the Wor…
6 Afghanistan 1954 CO2 Emissions (Metric Ton… Emissions 106 Rest of the Wor…
7 Afghanistan 1955 CO2 Emissions (Metric Ton… Emissions 154 Rest of the Wor…
8 Afghanistan 1956 CO2 Emissions (Metric Ton… Emissions 183 Rest of the Wor…
9 Afghanistan 1957 CO2 Emissions (Metric Ton… Emissions 293 Rest of the Wor…
10 Afghanistan 1958 CO2 Emissions (Metric Ton… Emissions 330 Rest of the Wor…
# … with 16,906 more rowsWe also need to sum the emissions across countries for each year. Here, we use the group_by() and summarize() function that we previously learned about.
data_long %>%
filter(Indicator == "Emissions") %>%
group_by(Year) %>%
summarize(Emissions = sum(Value))# A tibble: 264 x 2
Year Emissions
* <dbl> <dbl>
1 1751 9360
2 1752 9360
3 1753 9360
4 1754 9370
5 1755 9370
6 1756 10000
7 1757 10000
8 1758 10000
9 1759 10000
10 1760 10000
# … with 254 more rowsThen, we use the aes() argument of the ggplot() function to define that our x-axis will be the Year variable, the y-axis will be the emission Value variable, and that our data should be grouped or separated by the Country variable.
data_long %>%
filter(Indicator == "Emissions") %>%
group_by(Year) %>%
summarize(Emissions = sum(Value)) %>%
ggplot(aes(x = Year, y = Emissions))
Looks like we got a blank plot. What happened? We need to tell R what type of plot we want. To do that, we need to add another layer to define how we want the plot to look. We do so by using the + sign in between each command.
To tell R what type of plot we want, we need to add another layer to our ggplot object. To add a type of plot, we can use one of the many geom_* functions in ggplot2. For example, type geom into the RStudio console and you will see many options to scroll through.
Here, we will use the geom_line() function because we would like to create a line plot. We also use the size argument in geom_line() to control the size of the line.
data_long %>%
filter(Indicator == "Emissions") %>%
group_by(Year) %>%
summarize(Emissions = sum(Value)) %>%
ggplot(aes(x = Year, y = Emissions)) +
geom_line(size = 1.5)
Wow, the CO2 is really rising sharply!
Finally, let’s make this plot really nice by adding a few final touches. To change title, caption, and the axis labels, we can use the labs() function. Again, notice that a plus sign is used between each layer that we add to the plot. To make CO2 appear with a subscript we can use ~CO[2]~.
data_long %>%
filter(Indicator == "Emissions") %>%
group_by(Year) %>%
summarize(Emissions = sum(Value)) %>%
ggplot(aes(x = Year, y = Emissions)) +
geom_line(size = 1.5) +
labs(title = "World " ~CO[2]~ " Emissions per Year (1751-2014)",
caption = "Limited to reporting countries",
y = "Emissions (Metric Tonnes)")
Next, we use the theme() function to change the font size of the x-axis, y-axis, axis titles, and the caption as shown below. To know what to call each element of the plot in this function to change the size type ?theme() in the console.
You will see a very large list that includes other plot aspects like the background and the legend. This function can be used to modify your plot to your specifications.
data_long %>%
filter(Indicator == "Emissions") %>%
group_by(Year) %>%
summarize(Emissions = sum(Value)) %>%
ggplot(aes(x = Year, y = Emissions)) +
geom_line(size = 1.5) +
labs(title = "World " ~CO[2]~ " Emissions per Year (1751-2014)",
caption = "Limited to reporting countries",
y = "Emissions (Metric Tonnes)") +
theme_linedraw() +
theme(axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title.x = element_text(size = 12),
axis.title.y = element_text(size = 12),
plot.caption = element_text(size = 12),
plot.title = element_text(size = 16))
We can clearly see that global CO2 emissions have dramatically risen since 1900.
Notice, we used the function theme_linedraw() of ggplot2 to change the general appearance of the plot.
Useful tip: You can type theme_ in the RStudio console to see the various plot theme options available.
Great! We’ve created our first plot.
Before we leave this section, let’s save this theme so we do not have to keep typing the same code in future plots.
my_theme <-
theme_linedraw() +
theme(axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title.x = element_text(size = 12),
axis.title.y = element_text(size = 12),
plot.caption = element_text(size = 12),
plot.title = element_text(size = 16))In this way, we can just add another layer to our plot with the my_theme we created with the specifications of the sizes of the title and axis labels.
CO2_world <-
data_long %>%
filter(Indicator == "Emissions") %>%
group_by(Year) %>%
summarize(Emissions = sum(Value)) %>%
ggplot(aes(x = Year, y = Emissions)) +
geom_line(size = 1.5) +
labs(title = "World " ~CO[2]~ " Emissions per Year (1751-2014)",
caption = "Limited to reporting countries",
y = "Emissions (Metric Tonnes)") +
my_themeWe are also saving the plot to an object called CO2_world. To show the plot we simply type the name of the object:
Now let’s say we wanted to save this plot.
We could do so using the using the save() function to save this to a “plot” directory in our working directory as an RDA file and we can use the png() function to save a png for collaborators. We need to use dev.off() function to close the graphical device that we will use to create the png version of the plot so that we are ready to make another plot like this.
save(CO2_world, file =here::here("plots", "CO2_world.rda"))
png(here::here("plots", "CO2_world.png"))
CO2_world
dev.off()quartz_off_screen
2 One thing that would be nice to know is which countries are contributing the most or the least to CO2 emissions. Let’s continue to explore the data to investigate how CO2 emissions from individual countries have changed over time.
Next, we go back to using our data_long dataset. Here, we use the group argument in aes() which controls whether a line should be drawn
data_long %>%
filter(Indicator == "Emissions") %>%
ggplot(aes(x = Year, y = Value, group = Country)) +
geom_line() +
ylab("Emissions") +
my_theme
We can see that many countries show a dramatic increase in emissions over time with a handful of countries with particularly high levels.
Question Opportunity
group = Country argument?aes() arguments) can be changed?Since we have many overlapping lines, we will make our lines slightly transparent by using the alpha argument.
This takes values from 0 to 1, where 0 is completely transparent and 1 is completely opaque.
We also add our my_theme controlling the size of the title and axis labels.
CO2_countries <-
data_long %>%
filter(Indicator == "Emissions") %>%
ggplot(aes(x = Year, y = Value, group = Country)) +
geom_line(alpha = 0.4) +
labs(title = "Country" ~CO[2]~ "Emissions per Year (1751-2014)",
caption = "Limited to reporting countries",
y = "Emissions (Metric Tonnes)") +
my_theme
CO2_countries
Our plot is starting to look really good. One question you still might have is which country corresponds to which line? Which line indicates the emissions in the US?
We can add another “layer” on top of our first plot to add a blue line just for the US data. To do this we need to indicate what data we would like to plot, so we need to filter for just the US data and then we need to indicate that it will be colored by Country, even though in this case we only have one line to color. The default color would be a salmon pink color, but we would like blue. So we will use the scale_color_manual() function to manually choose the color that we want by using scale_colour_manual(values = c("blue")). Often you might use red to highlight a subset of the data, however, this can be difficult for viewers with certain types of colorblindness.
Notice how the color name needs to be in quotes and that the argument values = is used to specify what color values to use.
We can add this line to the plot in two ways.
Useful tip: Instead of retyping the original code to create the CO2_countries plot, we can just use the + operator:
CO2_countries +
geom_line(data = data_long %>%
filter(Indicator == "Emissions",
Country == "United States"),
aes(x = Year, y = Value, color = Country)) +
scale_colour_manual(values = c("blue"))
It looks like the US has long been the largest CO2 emission producing country until recently, when the US was surpassed by another country.
Let’s figure out who are the top 10 emission producing countries in 2014. Here, we filter the data for the year 2014, which was the final year of the data. Then, we can make a rank variable based on the Value variable for the amount of emissions produced.
There are many functions in the dplyr package for ranking values that are based on the SQL or specifically SQL rank functions. SQL is another programming language for managing large amounts of data. The difference in the rank functions mostly has to do with how to deal with ties in the data.
We will use dense_rank(), as we do not want gaps between ranks.
We want to do this in descending order because we want to rank by largest to smallest, so we will use the desc() function of the dplyr package. Then, we will arrange the output by rank using the arrange() function of the dplyr package.
top_10_count <-
data_long %>%
filter(Indicator == "Emissions", Year == 2014) %>%
mutate(rank = dense_rank(desc(Value))) %>%
filter(rank <= 10) %>%
arrange(rank)
top_10_count# A tibble: 10 x 7
Country Year Label Indicator Value Region rank
<fct> <dbl> <chr> <chr> <dbl> <chr> <int>
1 China 2014 CO2 Emissions (Metr… Emissions 1.03e7 Rest of the … 1
2 United Stat… 2014 CO2 Emissions (Metr… Emissions 5.25e6 United States 2
3 India 2014 CO2 Emissions (Metr… Emissions 2.24e6 Rest of the … 3
4 Russia 2014 CO2 Emissions (Metr… Emissions 1.71e6 Rest of the … 4
5 Japan 2014 CO2 Emissions (Metr… Emissions 1.21e6 Rest of the … 5
6 Germany 2014 CO2 Emissions (Metr… Emissions 7.20e5 Rest of the … 6
7 Iran 2014 CO2 Emissions (Metr… Emissions 6.49e5 Rest of the … 7
8 Saudi Arabia 2014 CO2 Emissions (Metr… Emissions 6.01e5 Rest of the … 8
9 South Korea 2014 CO2 Emissions (Metr… Emissions 5.87e5 Rest of the … 9
10 Canada 2014 CO2 Emissions (Metr… Emissions 5.37e5 Rest of the … 10We can see that China is now the top emission producing country.
Question Opportunity
What are the bottom 10 emission producing countries in 2014?
Let’s make a plot of just these top ten countries.
To do this, we need to filter the data to just these top countries by using the %in% operator to only keep countries in our Country variable that are also in the Country variable within top_10_count. We can use the pull() function also of the dplyr package to specifically grab just the Country data out of top_10_count.
Since we have 10 countries we will want to differentiate them by color.
To color our plot we will use the viridis color palette which is compatible with color-blindness by using the scale_color_viridis_d() function. This function is available by loading the ggplot2 package. There are a few variations for discrete values as _d, or binned continuous values as _b, or continuous values as _c. See here for more information.
Top10b <- data_long %>%
filter(Country %in% pull(top_10_count, Country)) %>%
filter(Indicator == "Emissions") %>%
filter(Year >= 1900) %>%
ggplot(aes(x = Year, y = Value, color = Country)) +
geom_line() +
scale_color_viridis_d() +
labs(title = "Top 10 Emissions-producing Countries in 2010 (1900-2014)",
subtitle = "Ordered by Emissions Produced in 2014",
y = "Emissions (Metric Tons)",
x = "Year") +
my_theme
Top10b
It’s still a bit difficult to tell which line corresponds to which country. So, let’s add a text label directly to the plot. .
One way to do this is to add text layer to our plot using the geom_text() function of the ggplot2 package. We need to first specify what variable or column we will use. However, we only want to pull the text labels for the top ten countries in the last year. To do this, we use the last() function of the dplyr package.
Then, we need to indicate that our text label will be based on the Country variable using the aes() aesthetics mapping argument. We will also get rid of our legend since we will not need it anymore, by using the theme() function of the ggplot2 package.
Top10b +
geom_text(data = data_long %>%
filter(Country %in% pull(top_10_count, Country)) %>%
filter(Indicator == "Emissions") %>%
filter(Year == last(Year)),
aes(label = Country)) +
theme(legend.position = "none")
Not bad, but some of the labels are overlapping and difficult to read. We can use the check_overlap = TRUE argument within the geom_text() function to remove overlapping variables. Also, we can expand the plot area horizontally so that the names are not cutoff by using scale_x_continuous(expand = c(0.2,0)). This takes a vector with two values. The first value indicates what percentage to expand the x axis in both directions. In our case we will expand by 15 percent. The second value indicates what absolute value to expand the limit of the x axis. In our case c(0.15,0) will achieve a similar result as c(0, 17), as the range of values from 1990 to 2014 is 114 years and 15% of this is 17.
Top10b +
geom_text(data = data_long %>%
filter(Country %in% pull(top_10_count, Country)) %>%
filter(Indicator == "Emissions") %>%
filter(Year == last(Year)),
aes(label = Country),
check_overlap = TRUE) +
scale_x_continuous(expand = c(0.15, 0)) +
theme(legend.position = "none")
This is easier to read now, but it also causes us to lose some of the labels. There are several alternative ways we can keep all of our labels and make them easier to read. The first package we will show is called directlabels.
The most simple option is to use the direct.label() function, which will automatically add labels at the end of the lines. However, it is a bit difficult to see some of our labels as they get automatically sized to fit the plot.
Alternatively this can be done within the ggplot2 framework by layering using the geom_dl() function.
Top10b +
scale_x_continuous(expand = c(0.3, 0)) +
geom_dl(aes(label = Country), method = list("last.bumpup")) +
theme(legend.position = "none")
This is more legible now. We have all 10 countries names listed and they are in order of the last data point and they are relatively close to the lines that they correspond to.
Another option is to use a different method in the directlables package. Here is a list of options.
For example, the "angled.boxes" method looks nice for some plots but does not work very well for our plot:
However the "last.polygons" method works quite well:
The second package is the ggrepel package, which is especially good for crowded labels that might overlap one another. It allows for more control than the directlabels package.
Specifically, we will use the geom_text_repel() function from the ggrepel package. Just like with geom_text(), first we need to specify what data we want to include. Then, we specify with the aes() argument that our label will be based on the Country variable and we again specify what variable to use for our x axis and y axis, so that we indicate where the labels should be plotted.
Top10b +
geom_text_repel(data = data_long %>%
filter(Country %in% pull(top_10_count, Country)) %>%
filter(Indicator == "Emissions") %>%
filter(Year == last(Year)),
aes(label = Country, x = Year, y = Value)) +
theme(legend.position = "none") +
scale_x_continuous(expand = c(0.3, 0))
You can see that this package creates segments that connect the label to the line.
There are many arguments to use to style your labels just the way that you want:
See the ggrepel vignette for more details.
Let’s play around with some of these options in the table above.
Top10b +
geom_text_repel(data = data_long %>%
filter(Country %in% pull(top_10_count, Country)) %>%
filter(Indicator == "Emissions") %>%
filter(Year == last(Year)),
aes(label = Country, x = Year, y = Value),
nudge_x = 10,
hjust = 1,
vjust = 1,
segment.size = 0.25,
force = 1) +
theme(legend.position = "none") +
scale_x_continuous(expand = c(0.3, 0)) +
scale_y_continuous(expand = c(0.3, 0))
Nice, that looks pretty good. For fun, let’s try showing our data in an entirely different way.
This time we will create a geom_tile plot.
To create this plot we will filter our data to include only the Countries included in the Country variable of the top_10_count.
Then, we will use the fct_reorder() function of the forcats package to order our countries based on the last emission value in 2014.
To use this function, the variable that is to be reordered is listed first. The variable that is being used to determine the order is listed second. Finally, a function to apply to the variable listed second is listed third. This function is used to determine the order. In this case, we want to determine the last value of the Value variable using the last() function (recall that this is also a function of the dplyr package). Then, the Country variable will be ordered by the last value of the Value variable.
To color our plot we will use the viridis color palette again but this time we will use the scale_fill_viridis_c() – recall that the _c indicates a continuous scale. See the scale_viridis reference for more information.
Top10t <-
data_long %>%
filter(Country %in% pull(top_10_count, Country)) %>%
filter(Indicator == "Emissions") %>%
filter(Year >= 1900) %>%
ggplot(aes(x = Year, y = fct_reorder(Country, Value, last))) +
geom_tile(aes(fill = log(Value))) +
scale_fill_viridis_c()Finally, let’s clean up the axes and the axes labels:
Top10t <- Top10t +
scale_x_continuous(breaks = seq(1900, 2014, by = 5),
labels = seq(1900, 2014, by = 5)) +
labs(title = "Top 10 " ~CO[2]~ "Emission-producing Countries",
subtitle = "Ordered by Emissions Produced in 2014",
fill = "Ln(CO2 Emissions (Metric Tonnes))") +
theme_classic() +
theme(axis.text.x = element_text(size = 12, angle = 90, color = "black"),
axis.text.y = element_text(size = 12, color = "black"),
axis.title = element_blank(),
plot.caption = element_text(size = 12),
plot.title = element_text(size = 16),
legend.position = "bottom")
Top10t
Now let’s say we wanted to save this plot.
We could do so using the using the save() function to save this to a “plot” directory in our working directory as an RDA file and we can use the png() function to save a png for collaborators. We need to use dev.off() function to close the graphical device that we will use to create the png version of the plot so that we are ready to make another plot like this.
save(Top10t, file =here::here("plots", "Top10t.rda"))
png(here::here("plots", "Top10t.png"))
Top10t
dev.off()quartz_off_screen
2 We see that Germany had very low emission rates at the end of World War II. We also see that the US has consistently had high emission rates since 1900, but the emission rates in China recently surpassed that of the US. The portions of the plot that are white indicate that there is no emission data for that country.
Question Opportunity
Think about what the pros and cons are of tile and line plots. In what situations would a tile plot be a better choice for effective scientific communication? In what situations would a line plot be a better choice?
Now, we will visualize all the variables in our dataset.
Here, we use the facet_wrap() function of the ggplot2 package, which plots multiple subplots simultaneously.
To use facet_wrap() with the option for a different y-axis scale for each subplot, we need to set the scales argument equal to "free_y". We can also indicate where we would like the label for the subplots to be located by using the strip.position argument.
ggplot(data_long, aes(x = Year, y = Value, group = Country)) +
geom_line(alpha = 0.2) +
geom_line(data = data_long %>%
filter(Country == "United States"),
aes(x = Year, y = Value, color = Country)) +
scale_colour_manual(values = c("blue")) +
labs(title = "Distribution of Indicators by Year and Value",
y = "Indicator Value") +
my_theme +
theme(strip.text = element_text(size = 16, face = "bold")) +
facet_wrap(Indicator ~ .,
scales = "free_y",
strip.position = "right",
ncol = 1)
Notice that we can change the size or style of the font for these labels using the strip.text = argument of the theme() function. We can also specify how many rows or columns we would like the subplots to be shown.
We can also facet by more than one variable (e.g. Indicator and Region to show the data from the US compared to other countries).
In this case we want the same y-axis to be used across the rows. We will use the facet_grid() function this time instead of facet_wrap() because of the way that the two facet variables are displayed. The facet_grid() function will as you might expect, create an output of plots that are displayed in a grid.
The syntax here is to put the name of the two variables on the left or right side of the ~ symbol, which tells you to facet by rows (left) or columns (right).
First, we will filter out the data about disasters and temperature as this is only for the US, by using the filter() function and ! indicates that we want only values of the Indicator variable not in the list containing "Disasters" and "Temperature".
data_long %>%
filter(!(Indicator %in% c("Disasters", "Temperature"))) %>%
ggplot(aes(x = Year, y = Value, group = Country)) +
geom_line() +
facet_grid(Indicator ~ Region, scales = "free_y") +
labs(title = "Distribution of Indicators by Year and Value",
y = "Indicator Value") +
my_theme +
theme(strip.text = element_text(size = 16, face = "bold"))
From these plots we can see that each type of data spans a different time span.
Question Opportunity
What happens when you create the same plot with facet_wrap()? Why might this be preferable in certain cases?
Question Opportunity
Calculate the total number of countries per year reporting C02 emissions, energy use and GDP. Plot this summary statistic (y-axis) across the years (x-axis). What do you see?
Hint : Use the tally() function in the dplyr package.
There are also some other common visualization techniques that are good for showing the difference between a set of observations and a mean value across time.
One of those is a line segment plot. For simplicity, let’s focus only on the data from the US. Let’s calculate the mean across all years for the CO2 emission and temperature from 1980 to 2010. Recall that our Indicator variable describes what kind of data we have (Emissions, Temperature, GDP, Energy, Disasters). We will calculate the mean for each Indicator set of data by first grouping by this variable and then calculating the mean of the values of the Value variable for each set of data. We will call this new variable Mean. Then, we calculate the difference between each observation and the mean and create a new variable for these values called Diff_from_mean. Again, all of this is performed for each group of data separately. Once we have created the new variables, we want to use the ungroup() function so that we no longer perform functions on subsets of the data based on the Indicator variable. Finally, we will also create a factor variable about the sign of the Diff_from_mean value to distinguish positive or negative changes. We will use this to color our plots.
data_long_us <-
data_long %>%
filter(Country == "United States", Year >= 1980, Year <= 2010) %>%
group_by(Indicator) %>%
mutate(Mean = mean(Value), Diff_from_mean = Value - Mean) %>%
ungroup() %>%
mutate(Diff_color = sign(Diff_from_mean)) %>%
mutate(Diff_color = as.factor(Diff_color))Rows: 155
Columns: 9
$ Country <fct> "United States", "United States", "United States", "Uni…
$ Year <dbl> 1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987, 1988, 1…
$ Label <chr> "CO2 Emissions (Metric Tons)", "CO2 Emissions (Metric T…
$ Indicator <chr> "Emissions", "Emissions", "Emissions", "Emissions", "Em…
$ Value <dbl> 4720000, 4540000, 4310000, 4340000, 4480000, 4490000, 4…
$ Region <chr> "United States", "United States", "United States", "Uni…
$ Mean <dbl> 5134194, 5134194, 5134194, 5134194, 5134194, 5134194, 5…
$ Diff_from_mean <dbl> -414193.548, -594193.548, -824193.548, -794193.548, -65…
$ Diff_color <fct> -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,…Next, we use the geom_segment() function to draw a straight line between points (x, y) and (xend, yend). In our case, this creates a plot that shows a bar for the difference between the observation and the mean across all the years.
data_long_us %>%
filter(Indicator %in% c("Emissions", "Temperature", "Disasters")) %>%
ggplot(aes(x = Year, y = Value)) +
geom_segment(aes(x = Year, y = Value, xend = Year,
yend = Mean, color = Diff_color),
size = 3.25) +
scale_color_manual(values = c("blue", "red")) +
geom_hline(aes(yintercept = Mean), linetype = 1, color = "black") +
facet_wrap(Indicator ~ ., scales = "free_y", ncol = 1) +
theme_classic() +
theme(axis.text.x = element_text(angle = 90),
axis.title = element_blank(),
legend.position = "none") +
labs(title = "US Disasters, Emissions, and Temperatures (1990-2010)",
subtitle = "Indicator Mean of 1990-2010 Represented by Solid Black Line")
We can see from this plot that overall there has been an increase in disasters, emissions and temperature in the most recent years.
Question Opportunity
What trends do you see in GDP and energy use across time?
Next, let’s zoom in on two of the variables: CO2 emissions and temperature. We use years between 1980 and 2014 as we have values for all of those years for these two variables.
We know that the datasets do not span the same amount of time. So let’s limit this plot to only the years where the data overlaps for both CO2 emissions and temperature.
We use the geom_point() function to create a scatter plot between the x and y variable defined in aes() where x is time and y is one of the two variables. We also add a line on top of the scatter plot that smooths the trend from the points. The smoother we used here is loess locally estimated scatterplot smoothing, a type of local polynomial regression fitting. This is a nonparametric regression that is also called a “moving regression” where subsets of points that are close to one another (hence the term local) are used in a least squares linear or nonlinear fit. Thus this results in a fit that may curve with the data.
CO2_temp_US_facet <-
data_long %>%
filter(Country == "United States", Year >= 1980, Year <= 2014,
Indicator %in% c("Emissions", "Temperature")) %>%
ggplot(aes(x = Year, y = Value)) +
geom_point() +
geom_smooth(method = "loess", se = FALSE) +
scale_x_continuous(breaks = seq(1980, 2014, by = 5),
labels = seq(1980, 2014, by = 5)) +
facet_wrap(Label ~ ., scales = "free_y", ncol = 1) +
theme_classic() +
theme(axis.text.x = element_text(size = 12, angle = 90, color = "black"),
axis.text.y = element_text(size = 12, color = "black"),
strip.text.x = element_text(size = 14),
axis.title = element_blank(),
plot.title = element_text(size = 16))
labs(title = "US Emissions and Temperatures (1980-2014)")$title
[1] "US Emissions and Temperatures (1980-2014)"
attr(,"class")
[1] "labels"
Note, we are showing a different theme here, namely the theme_classic() theme. We can see that there are similar patterns of CO2 emission levels and average annual temperatures.
We will save this plot now like so to our “plots” directory:
save(CO2_temp_US_facet, file =here::here("plots", "CO2_temp_US_facet.rda"))
png(here::here("plots", "CO2_temp_US_facet.png"))
CO2_temp_US_facet
dev.off()quartz_off_screen
2 Question Opportunity
geom_smooth().Next, instead of looking at the variables separately in faceted plots, let’s look at the relationship between CO2 emissions and other variables directly. Thus, it is useful to have each of these indicators as their own variable.
We can do this by using pivot_wider() to transform our long data table into a wide format.
wide_US <-
data_long %>%
filter(Country == "United States", Year >= 1980, Year <= 2014) %>%
select(-Label) %>%
pivot_wider(names_from = Indicator, values_from = Value)Let’s save this data as an rda file for future use and as a csv file, as this is often useful for collaborators. We will save this in a “wrangled” subdirectory of our “data” directory of our working directory.
save(wide_US, file = here::here("data", "wrangled", "wrangled_US_data.rda"))
readr::write_csv(wide_US, path = here::here("data","wrangled", "wrangled_US_data.csv"))Now we can specify which indicators we want to look at, so now we can specifically look at emissions and temperature.
CO2_temp_US <-
wide_US %>%
ggplot(aes(x = Emissions, y = Temperature)) +
geom_point() +
theme_classic() +
theme(axis.text.x = element_text(size = 12, color = "black"),
axis.text.y = element_text(size = 12, color = "black"),
axis.title = element_text(size = 14),
plot.title = element_text(size = 16)) +
labs(title = "US Emissions and Temperature (1980-2014)",
x = "Emissions (Metric Tonnes)",
y = "Temperature (Fahrenheit)")
CO2_temp_US
It might be helpful to add a trend line to this. We can do so by using the geom_smooth() function of the ggplot2 package.
If we want to look at a linear trend we need to specify the method using the method = lm argument. This adds a line to the data based on a linear model of the data using the lm function of the stats package. We will discuss the se = FALSE argument later.
We can just add this to the plot object that we just created to create a plot with this trend line.
Indeed, it does look like there is a positive, linear trend.
We will also save this plot:
save(CO2_temp_US, file =here::here("plots", "CO2_temp_US.rda"))
png(here::here("plots", "CO2_temp_US.png"))
CO2_temp_US
dev.off()quartz_off_screen
2 Question Opportunity
Now that we see that there might be a linear relationship between CO2 emissions and temperature, let’s learn about some statistical techniques to measure the strength of that relationship.
Suppose that we have a dataset called weight. Write some code to plot weight change by individual.
We want the plot to:
Note: first we need to use as.character() to convert “ID” from number to character.
If you are following along and stopped you could load the data you will need like so:
If you skipped the previous sections click here.
An RDA file (stands for R data) of the data (called wrangled_US_data.rda) can be found here or slightly more directly here. Download this file and then place it in your current working directory within a subdirectory called “wrangled” within a subdirectory called “data” to copy and paste our code. We used an RStudio project and the here package to navigate to the file more easily.
Click here to see more about creating new projects in RStudio.
You can create a project by going to the File menu of RStudio like so:
You can also do so by clicking the project button:
See here to learn more about using RStudio projects and here to learn more about the here package.
In this section, we are going to introduce some ways to better understand how two variables move together. We will focus on the CO2 emissions and temperature, but you will be encouraged to explore the relationship between CO2 emissions and the other variables.
We can always calculate the sample mean and variance for two variables.
# A tibble: 1 x 4
`mean(Emissions)` `mean(Temperature)` `sd(Emissions)` `sd(Temperature)`
<dbl> <dbl> <dbl> <dbl>
1 5142286. 52.9 450549. 0.891These are useful, but on their own they do not summarize whether or not there is a relationship between Emissions and Temperature (also these are on different scales entirely).
What else could we use? Next, we are going to learn about the correlation coefficient, which is a summary statistic that describes how two variables are related or move together.
We can use the correlation coefficient. Here, we are using this summary statistic to measure the strength of a linear relationship between two variables.
If we plot one variable on the x-axis and the other variable on the y-axis, we can see:
If the variables point upward in a very clear line, then there is a strong positive relationship. If the points do not really form a line, then there is a weak linear relationship or no linear relationship. There may however be a nonlinear relationship if the points create a different but defined shape.
See here for more information on nonlinear relationships.
If the points form a downward sloping line, then there is a negative relationship.
The numbers below each plot above are called correlation coefficients. They range from -1 to 1. A value of zero indicates that there is no correlation between the variables. While a value of 1 or -1 indicates perfect correlation, the closer the coefficient is to 1 or -1, the stronger the relationship. The sign of the coefficient indicates the direction of the relationship. If there is a negative relationship then the variables show opposing changes from each other - as one gets larger the other gets smaller. If the sign is positive, then the variables increase similarly.
Recall the plot we made above:
Let’s calculate the Pearson’s correlation coefficient called “Rho” \(\rho\) between CO2 emissions and temperature in the US. There are a few ways to calculate a correlation coefficient and this is one of the most common.
Formally, if we have a pair of observations \((x_1, y_1), \dots, (x_n,y_n)\), the correlation coefficient \(\rho\) between \(x\) and \(y\) is defined as
\[ \rho = \frac{1}{n-1} \sum_{i=1}^n \left( \frac{x_i-\mu_x}{\sigma_x} \right)\left( \frac{y_i-\mu_y}{\sigma_y} \right) \] where \(\mu_x, \mu_y\) are the means of \(x_1,\dots, x_n\) and \(y_1, \dots, y_n\), respectively, and \(\sigma_x, \sigma_y\) are the standard deviations.
Therefore, we can standardize the two variables and essentially average (the denominator is n-1) the standardized values to calculate the correlation coefficient rho.
Here we will manually perform the calculation. We will use the tally() function of the dplyr package to get the number of samples \(n\). In this case this is equivalent to the number of rows in the wide_US tibble. We need to then use the pull() function to specifically grab the value out of the tibble that is created from using this function. As you can see from using the base class() function that this is a tbl_df which is short for tibble data frame (the tidyverse version of a data frame) rather than just a number.
# A tibble: 1 x 1
n
<int>
1 35[1] "tbl_df" "tbl" "data.frame"When we check the class after using the pull() function we see that it is an integer.
[1] 35[1] "integer"We will also use the base scale() function to standardize the Emissions and Temperature values.
wide_US %>%
summarize(rho = (1/(pull(tally(wide_US), n) -1)) *(sum(scale(Emissions) * scale(Temperature)))) %>%
pull(rho)[1] 0.4711717Alternatively, you can use the cor() function in base R like so:
[1] 0.4711717If you want to learn more about why this is the calculation to determine the strength of the relationship between two variables, see this link.
Question Opportunity
There are different types of correlation coefficients. Look at the help file for the cor() function by typing ?cor() into the RStudio Console to learn more and then try a different method.
What are the differences?
Question Opportunity
To test if the association between a pair of variables is statistically significant, you can use the cor.test() of the stats package to calculate the correlation coefficient, as well as confidence intervals for correlation coefficients.
We can use the tidy() function of the broom package to make the output more usable for working with in R. The role of this function is to pull numeric values from outputs and create a data frame of the values.
Pearson's product-moment correlation
data: pull(wide_US, Emissions) and pull(wide_US, Temperature)
t = 3.0686, df = 33, p-value = 0.004277
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.1636155 0.6952523
sample estimates:
cor
0.4711717 # A tibble: 1 x 8
estimate statistic p.value parameter conf.low conf.high method alternative
<dbl> <dbl> <dbl> <int> <dbl> <dbl> <chr> <chr>
1 0.471 3.07 0.00428 33 0.164 0.695 Pearson's… two.sided We see that the correlation coefficient quantifying the strength of the linear relationship between C02 emissions and temperature is statistically significant.
Let’s briefly discuss the relationship between correlation and linear regression, which is further described in the Introduction to Data Science book.
We can use a regression line to predict a random variable \(Y\) given that we have gathered or observed some data about another variable \(X=x\). The regression line formally is defined as:
\[ \left( \frac{Y-\mu_Y}{\sigma_Y} \right) = \rho \left( \frac{x-\mu_X}{\sigma_X} \right) \] where \(\mu_X\) and \(\sigma_X\) (\(\mu_Y\) and \(\sigma_Y\)) are the mean and standard deviation of \(X\) (\(Y\)), and \(\rho\) is correlation between \(X\) and \(Y\). If \(x\) is larger than \(\mu_X\), then for every \(\sigma_X\), then \(Y\) will also increase \(\rho\) standard deviations above \(\mu_Y\).
Re-organizing the terms so that \(Y\) is on the left side and everything else is on the right side, we get:
\[ Y = \mu_Y + \rho \left( \frac{x-\mu_X}{\sigma_X} \right) \sigma_Y \] Thinking about some extreme examples:
To add regression lines to plots, we will need the above formula in the form:
\[ y= b + mx \mbox{ with slope } m = \rho \frac{\sigma_y}{\sigma_x} \mbox{ and intercept } b=\mu_y - m \mu_x \]
In our example, we can calculate the slope and intercept using the formula above and plot the line.
wide_US_summary <-
wide_US %>%
summarize(mu_x = mean(Emissions), sd_x = sd(Emissions),
mu_y = mean(Temperature), sd_y = sd(Temperature),
rho = cor(Emissions, Temperature),
slope = rho * sd_y / sd_x,
intercept = mu_y - rho * sd_y / sd_x * mu_x)
wide_US %>%
ggplot(aes(x = Emissions, y = Temperature)) +
geom_point() +
geom_abline(slope = wide_US_summary$slope,
intercept = wide_US_summary$intercept) +
my_theme
Note: In the plot above, we use the scale of the original variables (CO2 emission and temperature), but the formula above implies that standardization of the variables (i.e. subtracting the mean and dividing by the standard deviation) allows the regression line to have an intercept of 0 and slope equal to \(\rho\). A similar plot in standardized units is given below.
CO2_temp_US_sclaed <- wide_US %>%
ggplot(aes(x = scale(Emissions), y = scale(Temperature))) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
theme_linedraw() +
theme(axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14),
plot.caption = element_text(size = 12),
plot.title = element_text(size = 16))
labs(title = "US" ~ CO[2]~ "Emissions and Temperature (1980-2014)",
y = "Scaled Temperature (Fahrenheit)",
x = "Scaled Emissions (Metric Tonnes)")$y
[1] "Scaled Temperature (Fahrenheit)"
$x
[1] "Scaled Emissions (Metric Tonnes)"
$title
"US" ~ CO[2] ~ "Emissions and Temperature (1980-2014)"
attr(,"class")
[1] "labels"
Notice that we also use the geom_smooth(method = "lm") function that we previously used. Again, this adds a line corresponding to the slope and intercept from the lm() function from the stats R package.
Question Opportunity
What does se = FALSE mean? Try turning it to TRUE. What happens?
se stands for standard error. The gray shading shows the confidence interval of the smooth line.
While correlation is useful in many settings to understand how two variables move together, correlation is not always a useful summary. For example, here are ways it might not be useful:
As you can see from this plot, often data may show a similar pattern over time by random chance. See this website for more examples.
In this example and in our case study, data was collected over time. This type of data will often have what is called autocorrleation or serial correlation. This means that data points from one year to the next may be similar to one another or have some sort of internal structure related to time (such as seasonality). Let’s think about our CO2 emission and temperature data. You may be able to see how one year of CO2 emissions might be fairly similar to next year, because the number of sources (such as industrial factories and cars) in the US will change slightly from year to year, but they will be related to that of the previous years. Anytime we look at correlation between two variables that each have autocorrelation, this can result in a higher likelihood of correlation between these variables.
Indeed if we look at correlation between two random variables with autocorrelation we can see an inflation in the correlation rho values between the two variables, as compared to that of variables that do not have autocorrelation.
This is something to keep in mind when you evaluate how two things are related to one another over time. There are methods that have been developed to account for this in time series analysis (analyzing data over time). This does not mean that two variables (such as CO2 and temperature) may not be in fact correlated, it just means that we need to account for the autocorrelation within each variable, however this is beyond the scope of this case study.
The last thing we will do here is to create a plot that summarizes our major findings. We will use the plot_layout() function of the patchwork R package. The patchwork allows you to create a plot layout based on mathematical-like formulas. As you can see in this example we want the CO2_world and Top10t plot on top and we want another row with the CO2_temp_US_facet and the CO2_temp_US plot on the bottom. The plot_layout() function of the patchwork package then allows us to specify heights and widths for the plots.
We will also save the figure using the grDevices png() function. We can specify the name of our plot file and where we want it to be saved using the here() function of the here package. In this case in the img subdirectory of the directory containing our .Rproj file. We can also specify the size of the plot and the resolution using the res argument. The grDevices dev.off() function is necessary to close the graphics device.
This will include a few plots that we made previously in other sections. We saved these in the “plots” directory and will load these now for users who stopped and restarted or started at the data analysis section.
load(here::here("plots", "CO2_world.rda"))
load(here::here("plots", "Top10t.rda"))
load(here::here("plots", "CO2_temp_US_facet.rda"))png(here::here("img", "mainplot.png"), units = "in", width = 12, height = 10, res = 300)
(CO2_world | Top10t) / (CO2_temp_US_facet | CO2_temp_US) +
plot_layout(widths = c(1, 2),
heights = unit(c(4, 10), c('cm', 'cm')))
dev.off()
In this case study we evaluated CO2 emissions from as far back as 1751 for some countries to 2014. We discovered that global levels of CO2 emissions have dramatically increased over time. We also learned that some countries have been responsible for particularly high levels.
We also took a look at how CO2 emissions might relate to other factors, such as temperature, energy use, and natural disasters. We learned that we can summarize the relationship between two sets of data using correlation coefficients. We also learned that we can use regression to predict or describe how changes in one variable may influence changes in another variable. Importantly, we also learned that just because two variables show strong correlation or show an association, it does not necessarily indicate that they are causally related.
However, there is quite a bit of scientific evidence to indicate that in fact CO2 emissions trap heat and lead to increased global temperatures. Yet, it is important to realize that there are other factors involved in the relationship between US CO2 emissions and US annual average temperatures. For example there are CO2 emissions from other countries in the atmosphere, there are other greenhouse gases, there is already existing CO2 in the atmosphere that will continue to trap heat for many years, and finally there is heat trapped in the ocean due to previous emissions that will cause delayed changes in surface temperatures. However, it is vital that we work around the globe to reduce future greenhouse gas emissions to mitigate the increased temperatures that we will experience due to previous and existing CO2 emissions, so that the warming temperatures aren’t as extreme as they could be. Furthermore, we need to prepare for increased rates of natural disasters and how these may influence people around the world. Evidence suggests that impoverished people are the most affected by disasters. We need to be particularly mindful of this as we prepare for the future.
Ask students to create a plot with labels showing the countries with the lowest CO2 emission levels.
Ask students to plot CO2 emissions and other variables (e.g. energy use) on a scatter plot, calculate the Pearson’s correlation coefficient, and discuss results.
Tidyverse
RStudio cheatsheets Introduction to correlation Correlation coefficient
Correlation does not imply causation
Regression
Locally estimated scatterplot smoothing
Local polynomial regression
Autocorrleation
Time series
Methods to account for autocorrelation
US Environmental Protection Agency (EPA) Inventory of U.S. Greenhouse Gas Emissions and Sinks 2020 Report
National Climate Assessment Report
Greenhouse gases Climate change
Packages used in this case study:
| Package | Use in this case study |
|---|---|
here |
to easily load and save data |
readxl |
to import the excel file data |
readr |
to import the csv file data |
dplyr |
o view and wrangle the data, by modifying variables, renaming variables, selecting variables, creating variables, and arranging values within a variable |
magrittr |
to use and reassign data objects using the %<>%pipe operator |
stringr |
to select only the first 4 characters of date data |
purrr |
to apply a function on a list of tibbles (tibbles are the tidyverse version of a data frame) |
tidyr |
to drop rows with NA values from a tibble |
forcats |
to reorder the levels of a factor |
ggplot2 |
to make visualizations |
directlabels |
to add labels to plots easily |
ggrepel |
to add labels that don’t overlap to plots |
broom |
to make the output form statistical tests easier to work with |
patchwork |
to combine plots |
R version 4.0.3 (2020-10-10)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Mojave 10.14.6
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] patchwork_1.1.1 broom_0.7.4 ggrepel_0.9.1
[4] directlabels_2021.1.13 ggplot2_3.3.3 forcats_0.5.1
[7] tidyr_1.1.2 purrr_0.3.4 stringr_1.4.0
[10] magrittr_2.0.1 dplyr_1.0.4 readr_1.4.0
[13] readxl_1.3.1 knitr_1.31 here_1.0.1
loaded via a namespace (and not attached):
[1] tidyselect_1.1.0 xfun_0.21 bslib_0.2.4 lattice_0.20-41
[5] splines_4.0.3 colorspace_2.0-0 vctrs_0.3.6 generics_0.1.0
[9] viridisLite_0.3.0 htmltools_0.5.1.1 mgcv_1.8-33 yaml_2.2.1
[13] utf8_1.1.4 rlang_0.4.10 jquerylib_0.1.3 pillar_1.5.1
[17] glue_1.4.2 withr_2.4.1 DBI_1.1.1 lifecycle_1.0.0
[21] munsell_0.5.0 gtable_0.3.0 cellranger_1.1.0 evaluate_0.14
[25] labeling_0.4.2 ps_1.5.0 fansi_0.4.2 highr_0.8
[29] Rcpp_1.0.6 scales_1.1.1 backports_1.2.1 jsonlite_1.7.2
[33] farver_2.0.3 hms_1.0.0 digest_0.6.27 stringi_1.5.3
[37] rprojroot_2.0.2 grid_4.0.3 quadprog_1.5-8 cli_2.3.1
[41] tools_4.0.3 sass_0.3.1 tibble_3.1.0 crayon_1.4.1
[45] pkgconfig_2.0.3 Matrix_1.3-2 ellipsis_0.3.1 assertthat_0.2.1
[49] rmarkdown_2.6.6 rstudioapi_0.13 R6_2.5.0 nlme_3.1-151
[53] compiler_4.0.3 We would like to acknowledge Megan Latshaw for assisting in framing the major direction of the case study.
We would also like to acknowledge the Bloomberg American Health Initiative for funding this work.